Carolina Armella, Ágatha del Olmo y María Yu García
Published
October 18, 2024
1. Introducción
El propósito de este análisis es predecir el número esperado de reclamaciones que realizará una póliza a partir de la experiencia de siniestralidad histórica de una cartera de seguros. Este estudio se basa en la implementación de modelos lineales generalizados (GLM), que permiten realizar predicciones precisas manteniendo un alto grado de interpretabilidad.
Para este caso, utilizamos un conjunto de datos reales proporcionado por una compañía de seguros, que recoge información sobre 30,000 pólizas. En el presente análisis, para simplificar y reducir el tiempo de procesamiento, trabajamos con una muestra aleatoria de la base de datos original.
El análisis incluye una exploración detallada de los datos (Exploratory Data Analysis o EDA) y la implementación de modelos de aprendizaje automático para evaluar diferentes metodologías. Además, dado el contexto regulatorio de la industria de seguros, donde la transparencia y la explicabilidad son fundamentales, se priorizan los modelos que cumplen con estos requisitos, por eso no vamos a usar modelos de tipo caja negra.
2. Metodología
Para implementar un sistema automatizado de evaluación crediticia (credit scoring), es fundamental contar con modelos que permitan identificar con precisión los riesgos asociados a los préstamos. Además, la normativa bancaria en muchos países, incluido España, exige que las instituciones expliquen a los solicitantes los motivos para aceptar o rechazar una solicitud de crédito. Por ello, los modelos empleados deben ser transparentes. Herramientas más avanzadas como las redes neuronales o los modelos de máquinas de soporte vectorial (SVM) son consideradas “cajas negras” y no son fáciles (y a veces son imposibles) de interpretar y por lo tanto no cumplen con este requisito de transparencia.
2.1 Exploración de los datos
La base de datos corresponde a una extracción aleatoria de la BBDD de experiencia de siniestralidad de una compañía de seguros y consta de 19 variables capturando distintas características de las pólizas.
Reducción del dataset
Hemos decidido reducir la base de datos a 1000 observaciones en vez de las 27378 originales para agilizar la convergencia por problemas de rendimiento.
datos <-read.csv("D:\\3 curso\\TP\\Grupo 2.csv", sep=";")datos <- datos[1:1000,] # reducimos el training set para que tarde menosdim(datos)
En este caso, nuestra variable respuesta es numerosi, que indica el número de siniestros reportados por cada póliza. Esta variable es crucial porque el objetivo del modelo es predecir el número esperado de reclamaciones en función de las características de la póliza y del asegurado.
La variable nploiza es simplemente un identificador único de cada póliza, y no contiene información relevante o explicativa para la predicción de siniestros, por lo tanto la podemos eliminar.
Además, hay dos variables que nos dan información adicional. Podemos sustituir iniciop y finalp por la duración de la póliza en días o meses. Este predictor podría ser relevante teniendo en cuenta la exposición al riesgo. En cuanto al código postal, aunque en apariencia es solo una etiqueta numérica, contiene información geográfica valiosa, y para hacer más fácilmente interpretable este predictor extraemos la comunidad autónoma ya que no contamos con los metadatos de ccaa para saber a qué referencia cada número, por lo tanto esa variable la eliminamos.
2.2 Preprocesamiento de los datos
Eliminamos por las razones indicadas la variable npoliza.
# Descartamos el número de pólizadatos <-subset(datos, select =-c(npoliza))
Mostrar código
datos <- datos %>%mutate(region =substr(as.character(codigopostal), 1, 2)) %>%mutate(comunidad =case_when( region %in%c("01", "48", "20") ~"País Vasco", region %in%c("31") ~"Navarra", region %in%c("03", "12", "46") ~"Comunidad Valenciana", region %in%c("21", "41", "14", "29", "11", "18", "23", "04") ~"Andalucía", region %in%c("08", "17", "25", "43") ~"Cataluña", region %in%c("15", "27", "32", "36") ~"Galicia", region %in%c("33") ~"Asturias", region %in%c("39") ~"Cantabria", region %in%c("09", "42", "47", "34", "37", "05", "24", "40", "49") ~"Castilla y León", region %in%c("45", "13", "16", "19", "02") ~"Castilla-La Mancha", region %in%c("22", "44", "50") ~"Aragón", region %in%c("30") ~"Murcia", region %in%c("28") ~"Madrid", region %in%c("07") ~"Islas Baleares", region %in%c("35", "38") ~"Canarias", region %in%c("06", "10") ~"Extremadura", region %in%c("26") ~"La Rioja", region %in%c("51", "52") ~"Ceuta y Melilla" ))datos <- datos[, -19]
Comprobamos que no haya NULL en la variable creada ya que hemos contemplado todas las posibles combinaciones de códigos postales y sus comunidades autónomas según el INE en la fórmula.
sum(is.na(datos$comunidad))
[1] 260
Como vemos al crear la nueva variable comunidad autónoma, nos encontramos con muchos NA, por lo tanto, o la base de datos es completamente inventada o no sigue el sistema de códigos postales oficial, lo cual hace imposible su uso. Por lo tanto, eliminamos la variable creada.
datos <- datos[, -19]
Muchas de las variables no tienen la clase que les correspondería. Primero las variables iniciop y finalp las pasamos a formato fecha. También creamos la variable duración póliza con la diferencia entre estas dos, que nos indica cuánto ha durado la póliza en días.
A veces también es necesario convertir variables a tipo factor. Los factores se usan comúnmente para representar variables categóricas, es decir, variables que toman un número limitado de valores. En nuestro caso:
tipovehi: Representar diferentes tipos de vehículos, en concreto 7.
codigopostal: Se convierte en un factor, ya que representa códigos postales, que son categorías y no valores numéricos continuos.
prov: Similar, representa provincias que se pueden tratar como categorías.
ccaa: Representa la comunidad autónoma (en España), lo que también es una variable categórica.
Echamos un vistazo a las tablas de las variables no continuas para ver su frecuencia:
Creamos ahora una nueva variable llamada usovehi2 que clasifica la variable numérica usovehi en dos grupos, los que tienen un valor menor a 2000 que se pasan a llamar G1 y los que tienen un valor mayor, que se llamarán G2. Después eliminamos la variable original.
datos$usovehi2 <-ifelse(datos$usovehi <20000, "G1", "G2")# Descartamos uso de vehiculo, poca varianzadatos <-subset(datos, select =-c(usovehi))
Observamos la matriz de correlaciones, que nos indica el coeficiente de correlación entre todos los pares de variables numéricas de la base de datos.
Coeficiente de correlación
El coeficiente de correlación mide la intensidad y dirección de la relación lineal entre dos variables:
\(r=1\): Correlación positiva perfecta, a medida que una variable aumenta, la otra también aumenta proporcionalmente.
\(r=−1\): Correlación negativa perfecta, a medida que una variable aumenta, la otra disminuye proporcionalmente.
\(r=0\): No hay correlación lineal, las variables no tienen una relación lineal directamente.
numericas <- datos[sapply(datos, is.numeric)]# Calcular la matriz de correlacióncorrelaciones <-cor(numericas, use ="pairwise.complete.obs")# Crear un gráfico de correlacióncorrplot(correlaciones, method ="circle", type ="upper", tl.col ="black", tl.srt =45, col = colores)
Como podemos ver, lógicamente las variables están perfectamente correlacionadas de forma directa con ellas mismas. Pero podemos ver otras relaciones interesantes, para empezar hay una correlación bastante alta entre la variable objetivo numerosi y las variables cuantiasi y siniestro. Descartamos siniestro y cuantiasi como predictores.
cor(datos$siniestro, datos$numerosi)
[1] 0.7607234
cor(datos$cuantiasi, datos$numerosi)
[1] 0.3509457
datos <-subset(datos, select =-c(cuantiasi, siniestro))
Esto se debe a que la variable siniestro es un resultado de un evento ya ocurrido, y utilizarla para predecir el número de siniestros futuros es como usar la variable respuesta en el modelo de entrenamiento, lo que resulta en un modelo irrealista o sobreajustado por redundancias. Esto lo acabamos de ver reflejado en la matriz de correlaciones (\(r=0.76\)).
Además, al estar cuantiasi también altamente correlacionada con la variable de respuesta (\(r=0.35\)), usarla como predictor puede causar multicolinealidad por redundancia, lo que afecta negativamente la interpretación y precisión de los coeficientes del modelo.
El gráfico de correlaciones también nos ha mostrado una muy alta correlación entre potencia y valorvehi:
cor(datos$valorvehi, datos$potencia)
[1] 0.8207117
Esto nos lleva a un alto riesgo de multicolinealidad, que ocurre cuando dos o más variables predictoras están altamente correlacionadas entre sí, lo que puede causar problemas en los modelos estadísticos, especialmente en los modelos de regresión, ya que puede dificultar la interpretación de los coeficientes tal como hemos mencionado con cuantiasi. Por lo tanto, eliminamos la variable de menos interés interpretativo para el estudio, valorvehi, ya que la potencia es algo más objetivo que no depende de la economía.
datos <-subset(datos, select =-valorvehi)
Por último, el gráfico ha indicado una alta correlación entre edad y anticarnet.
cor(datos$anticarnet, datos$edad)
[1] 0.8174436
Si el estudio se enfocara en la relación entre estas dos variables para analizar el efecto de sacarse el carnet antes o después sería interesante mantener las dos, pero sabiendo que en la inmensa mayoría de casos las personas se sacan el carnet alrededor de los 18 años podríamos inferir de forma bastante acertada a partir de la edad la antigüedad del carnet. Por lo tanto, eliminamos la antigüedad de carnet.
datos <-subset(datos, select =-anticarnet)
Ahora podemos volver a ver el gráfico de correlaciones para comprobar que no quedan altas correlaciones:
numericas <- datos[sapply(datos, is.numeric)]correlaciones <-cor(numericas, use ="pairwise.complete.obs")corrplot(correlaciones, method ="circle", type ="upper", tl.col ="black", tl.srt =45, col = colores)
Al reducir el número de variables considerando las correlaciones, obtenemos una base de datos más manejable e interpretable, lo que facilita un análisis más eficiente y una convergencia más rápida en los modelos.
2.3 Separación de los datos en Training set y Test set
Partimos al azar base en dos: 66.7% entrenamiento (667 casos) y 33.33% test (333 casos).
Separación del dataset
Training set: Este conjunto se utiliza para entrenar el modelo. Contiene ejemplos históricos con características de préstamos y la información de si fueron pagados o no. El modelo analiza estos datos para identificar patrones que permitan clasificar a futuros solicitantes como buenos pagadores o posibles impagados.
Test set: Este segundo conjunto se emplea para evaluar la efectividad del modelo en datos que no ha visto antes. Sirve para comprobar si el modelo es capaz de reconocer patrones en casos nuevos y, por tanto, si sería útil para predecir con precisión en situaciones reales.
Para separar la base de datos se establece una semilla para garantizar que los resultados sean reproducibles. Luego, los datos se ordenan pseudo-aleatoriamente, se calcula el número de muestras para el conjunto de entrenamiento y de test, y se seleccionan los datos para ambos de la database original.
Después, se genera una tabla de frecuencias que muestra la distribución del número de reclamaciones en el conjunto de entrenamiento.
# Permutamos los datos al azarset.seed(216514)azar <-sample(1:nrow(datos))casos <-round(nrow(datos)*0.667)train <- datos[azar[1:casos],]test <- datos[azar[(casos+1):nrow(datos)],]Tabla.Training <-table(train$numerosi)Tabla.Test <-table(test$numerosi)Tabla.Training %>%kbl(col.names =c("Número Reclamaciones", "Frecuencia"),align =rep('c', 2)) %>%kable_styling()
En ambos casos la tabla sugiere que la mayoría de las pólizas en el conjunto de datos tienen 0 o 1 reclamación. Las otras categorías tienen menos frecuencia, lo que indica que las reclamaciones mayores son poco frecuentes en esta muestra, lo cual era de esperar, e interesa mantener a la empresa.
La técnica de separación usada es una de las más comunes y simples, pero puede no ser la mejor opción en todos los casos, especialmente si la distribución de las clases de los datos están desbalanceadas como acabamos de comprobar. Teniendo en cuenta esta característica de nuestra database, podemos usar el muestreo estratificado. Este método asegura que la distribución de las clases en los conjuntos de entrenamiento y test sea representativa de la distribución en el conjunto original. Es útil para evitar sesgos en la evaluación del modelo, especialmente cuando las clases menos frecuentes son importantes (que en nuestro caso lo son ya que son las que más afectan negativamente en la economía de la empresa).
3. Análisis Exploratorio de Datos
3.1 Variable respuesta
Visualizamos la frecuencia de número de reclamaciones a lo largo del tiempo.
Mostrar código
barplot(table(train$numerosi),col = colores, border ="white", main ="Frecuencia de Número de Reclamaciones", xlab ="Número de Reclamaciones", ylab ="Frecuencia", las =1, cex.names =0.9, cex.axis =0.8, cex.main =1.2)
Una de las propiedades más importantes de la distribución de Poisson es que su media es igual a su varianza. Ambos están representados por el parámetro λ (lambda). Esta propiedad es única y ayuda a identificar si un conjunto de datos sigue una distribución de Poisson. Esta igualdad implica que, a medida que aumenta el número esperado de sucesos, también lo hace la variabilidad en el número real de sucesos.
Por lo tanto, podemos encontraar sobredispersión si la varianza de la variable de respuesta es significativamente mayor que su media, lo cual es incompatible con los supuestos del modelo Poisson como acabamos de ver.
mean(train$numerosi)
[1] 0.4482759
var(train$numerosi)
[1] 0.9383867
La media es 0.45, mientras que la varianza es 0.94. En un modelo Poisson ideal, la media y la varianza deberían ser iguales.
La varianza es notablemente mayor que la media, lo que indica sobredispersión. Esto significa que los datos probablemente no sean bien modelados por un modelo Poisson estándar, ya que este asume que la media y la varianza son iguales. Esta igualdad implica que, a medida que aumenta el número esperado de sucesos, también lo hace la variabilidad en el número real de sucesos.
El exceso de ceros ocurre cuando la proporción de ceros en la variable de respuesta es mucho mayor de lo esperado bajo una distribución Poisson.
table(train$numerosi)[1] /nrow(train)
0
0.7331334
Más del 73% de las observaciones tienen un valor de cero en número de reclamaciones, lo cual esperaríamos porque indica que muchos pólizas no tienen reclamaciones, pero en términos técnicos nos puede llevar a complicaciones. Este exceso de ceros no es consistente con un modelo Poisson estándar, que generaría ceros con menor frecuencia dado el valor esperado de 0.39.
Concluimos que el modelo Poisson estándar no es adecuado. Necesitarás un modelo que pueda manejar sobredispersión y exceso de ceros como por ejemplo los modelos de inflación cero. Estos utilizan una técnica estadística adaptada a los datos de recuento que tienen un exceso de resultados cero, lo cual es ideal para nuestro caso.
Ahora, representamps el número de reclamaciones ajustadas según un modelo de Poisson, el número de reclamaciones que el modelo de Poisson predice que deberíamos tener para cada número de reclamaciones frente al número de reclamaciones observadas en los datos reales.
Mostrar código
predict0 <-dpois(0:10, mean(train$numerosi)) *nrow(train)observado <-table(factor(train$numerosi, levels =0:10))barplot(rbind(observado, predict0), las =1,beside =TRUE, names.arg =c(0:9, ">=10"),col =c(colores[1], colores[5]),border ="white", main ="Distribución Poisson", xlab ="Número de Reclamaciones", ylab ="Frecuencia", las =1, cex.names =0.9, cex.axis =0.8, cex.main =1.2)
Como podemos ver, el modelo Poisson predice que deberíamos tener menos pólizas con 0 reclamaciones que las observadas, pero más de 1 reclamación, lo cual implicaría una variabilidad más constante con la media.
Comparamos con una Binomial Negativa:
Mostrar código
fit.bn <-goodfit(train$numerosi, type ="nbinomial", method ="MinChisq")barplot(rbind(observado, fit.bn$fitted[1:11]), las =1,beside =TRUE, names.arg =c(0:9, ">=10"),col=c(colores[1], colores[5]),border ="white", main ="Distribución Binomial Negativa", xlab ="Número de Reclamaciones", ylab ="Frecuencia", las =1, cex.names =0.9, cex.axis =0.8, cex.main =1.2)
Claramente la distribución Binomial Negativa se ajusta de manera más adecuada a nuestra base de datos, ya que las frecuencias ajustadas se alinean mucho mejor con las observadas, reflejando de manera más precisa la estructura de los datos.
3.2 Relaciones variable respuesta-predictores
3.2.1 Tiempo de exposición al riesgo
El tiempo de exposición al riesgo se debe tratar de la siguiente manera según el contexto:
Como un offset: es el enfoque más común en modelos de conteos (Poisson, Binomial Negativa).
Como predictor: este enfoque es menos común, pero útil si creemos que el efecto del tiempo varía entre los casos.
Como peso: si queremos dar más importancia a los casos con más tiempo de exposición al riesgo.
Por lo tanto, para modelos de eventos o conteos, el tratamiento más adecuado es generalmente como un offset.
# Tiempo de exposición (Predictor, Peso u Offset)summary(train$duracion_poliza)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0 79.0 179.0 178.4 273.0 366.0
Como el tiempo de exposición afecta la tasa de siniestros, modelarlo como un offset nos permitiría ajustar por esta variable sin estimar un coeficiente independiente para ella. Al tratarlo como offset, estaríamos ajustando el modelo para que la tasa de eventos sea proporcional al tiempo de exposición, pero sin incluirla como un predictor adicional cuyo coeficiente se estimaría.
Es posible que algunos contratos tengan una duración extremadamente corta, de menos de un día, como un seguro de alquiler que se contrata y se extingue el mismo día. En este caso, la “exposición al riesgo” podría ser efectivamente cero, ya que no hay un período suficiente de tiempo para que ocurra un evento asegurado (como un accidente o siniestro).
sum(datos$duracion_poliza==0)
[1] 24
sum(train$duracion_poliza==0)
[1] 13
sum(test$duracion_poliza==0)
[1] 11
En nuestro caso, el conjunto de entrenamiento contiene 13 cero días al azar de los 24 registros con ceros originalmente. De la misma forma, en el conjunto de prueba hay 11 casos con duración de póliza igual a cero.
En ambos conjuntos hay registros con duración de póliza de 0 días, lo cual podría ser problemático, así que los convertimos en 1, y lo repetimos para la posibilidad de que haya valores erróneos menores a 0 que se deberían a un error de imputación de los datos originales.
Para entender mejor la distribución de la duración de las pólizas en nuestro conjunto de datos, hemos realizado un gráfico de densidad. Este gráfico nos permite observar cómo se distribuyen los valores de duración de las pólizas, destacando posibles tendencias o patrones en la variable.
Mostrar código
ggplot(train, aes(x = duracion_poliza)) +geom_density(fill = colores[1], color = colores[5], alpha =0.6) +theme_minimal() +labs(title ="Distribución de la Densidad de la Duración de la Póliza", x ="Duración de la Póliza (días)", y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm"))
La distribución de la duración de las pólizas muestra un comportamiento relativamente uniforme a lo largo del eje x, con el máximo alrededor de los 175 días. Los valores mínimos se observan al inicio y hacia el final del rango, indicando una menor frecuencia de pólizas con duraciones extremadamente cortas o largas.
Mostrar código
set.seed(1)uniforme <-runif(n =nrow(train), min =min(train$duracion_poliza), max =max(train$duracion_poliza))data_comparacion <-data.frame(duracion_poliza =c(train$duracion_poliza, uniforme),tipo =rep(c("Real", "Uniforme"), each =nrow(train)))ggplot(data_comparacion, aes(x = duracion_poliza, fill = tipo, color = tipo)) +geom_density(alpha =0.6) +theme_minimal() +labs(title ="Comparación de Distribución Duración Póliza y Uniforme", x ="Duración de la Póliza (días)", y ="Densidad") +scale_fill_manual(values =c(colores[1], colores2[1])) +scale_color_manual(values =c(colores[5], colores2[5])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),legend.title =element_blank(),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) +guides(fill =guide_legend(override.aes =list(alpha =0.6)))
Simulando una distribución uniforme aleatoria (con la semilla 1) podemos ver claramente que son en efecto bastante similares. Esto implica que todas las duraciones tienen una probabilidad similar de ocurrir, lo cual es consistente con una distribución uniforme, donde cada valor dentro del rango tiene la misma probabilidad.
Además, puede ser interesante estudiar el límite máximo de la variable:
max(train$duracion_poliza)
[1] 366
Que el valor máximo de la variable sea 366 días sugiere que el contrato abarca un máximo de un año calendario completo en un año bisiesto, aunque no es lo más común.
cor(train$duracion_poliza, train$numerosi)
[1] 0.2735884
Para interpretar los valores de correlación vamos a seguir estos intervalos:
Intervalos de interpretación del coeficiente de correlación:
\(|r|=1\): Correlación perfecta.
\(|r|\geqslant0.8\): Correlación alta.
\(|r|\geqslant0.5\): Correlación media
\(|r|\geqslant0.3\): Correlación débil
\(|r|>0\): Correlación muy débil
\(r=0\): No hay correlación.
La correlación entre la duración de la póliza y el número de siniestros es muy débil pero positiva (\(r=0.27\)). Esto sugiere que, en general, a medida que aumenta la duración de la póliza, también tiende a incrementarse el número de siniestros, aunque la relación no es especialmente fuerte. Una póliza más larga brinda más tiempo durante el cual puede ocurrir un siniestro. Por tanto, el aumento en el número de siniestros es razonable desde una perspectiva temporal.
De todas formas, aunque la relación existe, no es determinante. Otros factores además de la duración probablemente influyen significativamente en el número de siniestros, como el tipo de póliza, la antigüedad del vehículo o las características del asegurado. Podría ser interesante estudiar la media de siniestros para cada cantidad de días de duración.
Mostrar código
df_plot <-data.frame(duracion_poliza =sort(unique(train$duracion_poliza)),promedio_siniestros =tapply(train$numerosi, train$duracion_poliza, mean))ggplot(df_plot, aes(x = duracion_poliza, y = promedio_siniestros)) +geom_line(color = colores[1], size =1.2) +geom_point(color = colores[5], size =3, shape =19) +theme_minimal() +labs(title ="Promedio de Siniestros por Duración de la Póliza",x ="Duración de la Póliza (días)", y ="Promedio de Siniestros") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"), axis.title.x =element_text(size =12),axis.title.y =element_text(size =12),axis.text =element_text(size =10), panel.grid.major =element_line(color ="gray", size =0.5), panel.grid.minor =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
Cada punto en el gráfico muestra el promedio de siniestros para una duración específica de la póliza. Parece que hay más variabilidad en el promedio de siniestros a medida que aumenta la duración de la póliza. Esto significa que, conforme pasa el tiempo, es más común que las pólizas tengan diferentes números de reclamaciones. No se observa una hay una tendencia muy clara en los puntos del gráfico, lo que sugiere que el número promedio de siniestros no tiene un patrón fácil de predecir solo con la duración de la póliza.
3.2.2 Potencia del vehículo
summary(train$potencia)
Min. 1st Qu. Median Mean 3rd Qu. Max.
35.00 60.00 80.00 86.75 103.00 272.00
Mostrar código
ggplot(train, aes(x = potencia)) +geom_density(fill = colores[1], color = colores[5], alpha =0.6) +theme_minimal() +labs(title ="Distribución de la Densidad de la Potencia", x ="Potencia", y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm"))
Esta variable parece seguir una distribución similar a la log-normal. Lo comprobamos:
Mostrar código
set.seed(1)log_normal <-rlnorm(n =nrow(train), meanlog =mean(log(train$potencia +1)), sdlog =sd(log(train$potencia +1)))data_comparacion_potencia <-data.frame(potencia =c(train$potencia, log_normal),tipo =rep(c("Real", "Log-Normal"), each =nrow(train)))ggplot() +geom_density(data =subset(data_comparacion_potencia, tipo =="Real"), aes(x = potencia, fill ="Real"), alpha =0.6, color = colores[5]) +geom_density(data =subset(data_comparacion_potencia, tipo =="Log-Normal"), aes(x = potencia, fill ="Log-Normal"), alpha =0.4, color = colores2[5]) +theme_minimal() +labs(title ="Comparación de Distribución Potencia y Log-Normal", x ="Potencia", y ="Densidad") +scale_fill_manual(values =c("Real"= colores[1], "Log-Normal"= colores2[1])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),legend.title =element_blank(),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) +guides(fill =guide_legend(override.aes =list(alpha =0.6)))
La distribución real de la potencia y la simulada con una log-normal son muy similares, lo que sugiere que la potencia de los vehículos en el conjunto de datos podría aproximarse razonablemente a este tipo de distribución. Esto puede darse porque en muchas situaciones reales, las variables relacionadas con características físicas o de rendimiento (como la potencia del motor) tienden a distribuirse de forma asimétrica, con una mayoría de valores en un rango medio y una cola hacia valores más altos que en este caso serían los coches de carreras ode lujo. Este comportamiento es característico de una distribución log-normal.
Ahora estudiamoscómo varía el número promedio de siniestros (variable respuesta) en función de la potencia del vehículo:
plot(sort(unique(train$potencia)), medias, pch =20, col = colores[3], xlab ="Potencia del vehículo", ylab ="Número promedio de siniestros", main ="Relación entre Potencia y Número Promedio de Siniestros", cex =1.2, # Tamaño de los puntoscex.axis =1, # Tamaño de los ejescex.lab =1.1, # Tamaño de las etiquetascex.main =1.3, # Tamaño del títulolwd =2, # Grosor de la línea de los puntoslty =1, # Tipo de línea (sin línea, solo puntos)xlim =c(min(sort(unique(train$potencia))) -10, max(sort(unique(train$potencia))) +10))
Aparentemente, la mayoría de los puntos están agrupados en los valores más bajos de potencia y número promedio de siniestros, lo que sugiere que los vehículos con menor potencia tienden a tener menos siniestros. A parte, hay algunos puntos fuera de la masa principal de puntos, con valores más altos tanto en potencia como en el número promedio de siniestros. Esto indica que, aunque no es la norma, algunos vehículos con alta potencia tienen un número significativo de siniestros.
En nuestra base de datos no parece haber una tendencia específica, ya que los puntos se distribuyen al rededor de 0. Esto podría sugerir que los vehículos de menor potencia son más seguros o que son conducidos de manera más conservadora.
cor(train$potencia, train$numerosi)
[1] 0.08484551
Volviendo a la interpretación de cada rango de correlación, un valor de \(r=0.08\) indica una relación lineal directa muy débil, lo cual podíamos esperar del gráfico anterior ya que no se observa ninguna relación clara.
La discretización de una variable continua, como en este caso la potencia de un vehículo, implica dividirla en grupos o categorías. Este proceso puede ser útil para facilitar la interpretación de la relación entre esa variable y otras, como el número de siniestros. Se pueden usar los cuantiles de la variable para dividir los datos en grupos más equitativos en términos de cantidad de observaciones.
# Discretizar la potencia usando cuantilesqnt <-quantile(train$potencia, seq(0, 1, .2))qnt[1] <- qnt[1] -0.01# Ajustar el primer cuantílqnt[5] <- qnt[5] +0.01# Ajustar el quinto cuantíltrain$g.potencia <-cut(train$potencia, breaks = qnt)train$g.potencia <-factor(train$g.potencia, levels =levels(train$g.potencia))medias <-tapply(train$numerosi, train$g.potencia, mean)
Mostrar código
ggplot(data =data.frame(x =1:length(medias), y = medias), aes(x = x, y = y)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5) +scale_x_continuous(breaks =1:length(medias), labels =levels(train$g.potencia)) +labs(title ="Promedio de Siniestros por Grupo de Potencia", x ="Rango de Potencia", y ="Promedio de Siniestros") +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
Como podemos ver, los vehículos en el rango de potencia (90-112) parecen estar asociados con un mayor número de siniestros. Esto podría deberse a que estos vehículos tienen suficiente potencia para ser utilizados de manera más agresiva o en condiciones que aumentan el riesgo de siniestros.
Los vehículos en el rango de potencia (35-90) y (112-272) tienen menos siniestros en promedio, lo cual podría sugerir que estos vehículos son más seguros o que son conducidos de manera más conservadora, de hecho probablemente estos vehículos son tractores, patinetes y otros similares.
La variabilidad en otros rangos muestra que hay diferencias aparentemente significativas en el número promedio de siniestros según la potencia del vehículo.
La discretización de la variable ha sido muy interesante, podemos estar interesados en hacerlo de nuevo, por eso hemos creado una función para hacerlo de forma más automática:
Esta función nos ayudará a realizar los cálculos de discretización de forma mucho más rápida.
3.2.3 Antigüedad del vehículo
Estudiamos ahora la variable de antigüedad del vehículo.
Mostrar código
ggplot(train, aes(x = antivehi)) +geom_density(bw =0.5, fill = colores[3], color = colores[5], alpha =0.6) +labs(title ="Distribución de la Antigüedad del Vehículo",x ="Antigüedad del Vehículo (años)",y ="Densidad" ) +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12),axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
El gráfico no parece ajustarse a una distribución específica, pero sí muestra una clara tendencia decreciente. Observamos que hay muchos vehículos relativamente nuevos, con menos de 5 años de antigüedad, seguido de una disminución significativa a partir de entre los 11 y 12 años.
Para realizar la comparación, identificamos que la distribución más similar a nuestros datos es la chi-cuadrado, aunque no parece tan similar como en los casos anteriores. Esta distribución se caracteriza por su asimetría, comenzando con valores altos que disminuyen progresivamente. Es adecuada en este caso porque refleja una tendencia inicial elevada que va decreciendo, algo similar a lo observado en nuestros datos de antigüedad.
Mostrar código
scale_chi <-mean(train$antivehi) /2set.seed(1)chi_data <-rchisq(n =length(train$antivehi), df =2) * scale_chidata_comparacion_antivehi <-data.frame(antiguedad =c(train$antivehi, chi_data),tipo =c(rep("Real", length(train$antivehi)), rep("Chi-Cuadrado", length(chi_data))))ggplot() +geom_density(data =subset(data_comparacion_antivehi, tipo =="Real"), aes(x = antiguedad, fill ="Real"), alpha =0.6, color = colores[5]) +geom_density(data =subset(data_comparacion_antivehi, tipo =="Chi-Cuadrado"), aes(x = antiguedad, fill ="Chi-Cuadrado"), alpha =0.4, color = colores2[5]) +theme_minimal() +labs(title ="Comparación de Distribución Antigüedad del Vehículo y Chi-Cuadrado", x ="Antigüedad del Vehículo (años)", y ="Densidad") +scale_fill_manual(values =c("Real"= colores[1], "Chi-Cuadrado"= colores2[1])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),legend.title =element_blank(),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) +guides(fill =guide_legend(override.aes =list(alpha =0.6)))
Como podemos observar, ambas distribuciones son parecidas, salvo por un repunte en los datos reales alrededor de los 10 años de antigüedad. Este comportamiento podría estar relacionado con un mayor número de vehículos que alcanzan esa antigüedad debido a características del mercado o patrones de uso. A pesar de esta diferencia, la similitud se recupera al volver a descender posteriormente.
Mostrar código
antivehi_data <-data.frame(antivehi =sort(unique(train$antivehi)),mean_numerosi =tapply(train$numerosi, train$antivehi, mean))ggplot(antivehi_data, aes(x = antivehi, y = mean_numerosi)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5) +theme_minimal() +labs(title ="Relación entre Antigüedad del Vehículo y Número de Siniestros",x ="Antigüedad del Vehículo (años)",y ="Promedio de Número de Siniestros" ) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
Los vehículos más nuevos (0 años) tienen el promedio de siniestros más alto, más de 0.87. Esto podría deberse a la falta de experiencia del conductor con el vehículo nuevo. A medida que los vehículos son ligeramente más antiguos (de 1 a 5 años), el promedio de siniestros disminuye, alcanzando un mínimo local a los 7 años. Después, el promedio de siniestros muestra algunas fluctuaciones, pero en general vemos que decrece con algunos incrementos ocasionales.
Tras ver el gráfico, puede ser interesante ver la correlación entre la variable de respuesta, número de siniestros, y la antigüedad del vehículo.
cor(train$antivehi, train$numerosi)
[1] -0.1807625
Hay una muy débil correlación directa (\(r=-0.18\)), es decir, existe una pequeña tendencia a que, a medida que aumenta la antigüedad del vehículo, disminuye el número de siniestros en promedio. El hecho de que sea tan débil implica que la antigüedad del vehículo no es un predictor particularmente fuerte del número de siniestros.
3.2.4 Tipo de vehículo
La variable tipo de vehículo es de tipo factor, por lo tanto es preciso otro tipo de análisis.
summary(train$tipovehi)
100 120 150 200 250 300 310
581 8 32 31 2 3 10
Estudiamos su estructura con un gráfico:
Mostrar código
ggplot(train, aes(x = tipovehi, y = numerosi)) +stat_summary(fun = mean, geom ="bar", fill = colores[1], color = colores[5], width =0.6) +theme_minimal() +labs(title ="Promedio de Siniestros por Tipo de Vehículo", x ="Tipo de Vehículo", y ="Promedio de Siniestros") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),axis.text.x =element_text(angle =45, hjust =1), panel.border =element_blank(),plot.margin =margin(2, 1, 1, 1, "cm"))
El gráfico muestra el promedio de siniestros (reclamaciones) para diferentes tipos de vehículos. Los vehículos del tipo 120 y 310 tienen los promedios de siniestros más altos (aproximadamente 1 y 0.65 respectivamente). Esto sugiere que estos tipos de vehículos tienen más probabilidades de estar involucrados en siniestros. Los vehículos del tipo 250 y 300 no reportan a penas siniestros y los del resto de tipos tienen un promedio bajo (menos de 0.5). Esto indica que estos vehículos son menos propensos a involucrarse en siniestros.
Hay diferencias significativas en el promedio de siniestros entre los diferentes tipos de vehículos, lo cual es importante para los estudios de evaluación de riesgos en seguros. Estas diferencias hacen que la distribución no se asemeje a ninguna distribución común.
Ahora, al ser la variable de tipo factor, comprobamos si hay diferencias significativas en el número de siniestros entre los distintos tipos de vehículos. Usando ANOVA sabemos si tipo de vehículo tiene un impacto importante en el número de siniestros o si es más probable que las diferencias que vemos sean solo por casualidad.
Analysis of Variance Table
Response: train$numerosi
Df Sum Sq Mean Sq F value Pr(>F)
train$tipovehi 6 4.53 0.75521 0.8034 0.5675
Residuals 660 620.43 0.94005
Umbrales comunes de significancia
El coeficiente de correlación mide la intensidad y dirección de la relación lineal entre dos variables:
\(p-value=0.05\): Este es el umbral más comúnmente utilizado.
\(p-value=0.01\): Un umbral más estricto.
\(p-value=0.1\): A veces se usa en investigaciones exploratorias, pero no se considera concluyente.
El análisis de varianza muestra que el tipo de vehículo no tiene un impacto significativo en el número de siniestros. El valor F es bajo y el valor p es alto (0.5675), no entra en ninguno de los casos mencionados, lo que sugiere que las diferencias en el número de siniestros entre los diferentes tipos de vehículos podrían ser atribuibles al azar y no son estadísticamente significativas. En otras palabras, el tipo de vehículo no parece ser un factor importante para predecir el número de siniestros en este caso.
Aquí vemos los intervalos de confianza al 95% para los coeficientes estimados del modelo lineal que hemos creado. Específicamente, estos números nos dicen el rango dentro del cual creemos que se encuentra el valor real del coeficiente, con un 95% de certeza.
Para todos los tipos de vehículo mencionados, los intervalos de confianza incluyen el cero en su rango, lo que indica que, en efecto, no podemos afirmar que exista una relación significativa entre estos tipos de vehículos y el número de siniestros. En otras palabras, no podemos concluir que el tipo de vehículo tenga un impacto claro sobre el número de siniestros en el modelo.
3.2.5 Plazas de vehículo
Estudiamos ahora la variable plazas de vehículo.
table(train$plazasvehi)
5 6
662 5
Como vemos, es una variable muy desbalanceada, aunque en términos generales, plazasvehi tiene dos categorías (5 y 6 plazas), la gran mayoría de los vehículos tiene 5 plazas, mientras que solo una pequeña proporción tiene 6 plazas.
tapply(train$numerosi, train$plazasvehi, mean)
5 6
0.4471299 0.6000000
Esta funciónnos permite ver cómo varía el promedio de siniestros (numerosi) según el número de plazas del vehículo (plazasvehi), lo cual elimina el efecto del desbalance. El resultado nos indica que para los vehículos con 5 plazas, el promedio de los siniestros (train$numerosi) es 0.4471, y para los vehículos con 6 plazas, el promedio de los siniestros es 0.6.
Los vehículos con 6 plazas suelen ser más grandes, lo que podría implicar que son utilizados de manera diferente o para viajes más largos, lo que podría aumentar las probabilidades de que ocurran siniestros. En cambio, los vehículos más pequeños de 5 plazas a menudo son conducidos de manera diferente o menos veces, lo que podría reducir la exposición al riesgo de siniestros.
3.2.6 Edad
Ahora estudiamos la variable de edad. Realizamos de nuevo el gráfico de la distribución de la variable:
ggplot(train, aes(x = edad)) +geom_density(fill = colores[3], color = colores[5], alpha =0.6, bw =0.5) +theme_minimal() +labs(title ="Distribución de la Edad", x ="Edad", y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
La distribución de la edad sigue una forma que se asemeja bastante a una distribución normal, lo cual es común en los conjuntos de datos con esta variable. Las vemos constrastadas en el siguiente gráfico:
La distribución observada en los datos reales sigue una forma similar a una distribución normal, pero presenta algunas irregularidades. El gráfico muestra un pico notable cerca de los 50 años, lo que indica que la mayoría de los casos se concentran en este rango de edad. Además, se observa un descenso en la densidad justo después de los 50 años, lo que sugiere una disminución en la frecuencia de los datos en ese grupo, para el cual no queda muy clara la interpretación.
Estudiamos ahora su relación con la variable respuesta:
cor(train$edad, train$numerosi)
[1] -0.008924572
Para empezar, el coeficiente de correlación ($r=-0.008) nos indica una correlación muy débil negativa. A medida que la edad del conductor aumenta, el número de siniestros disminuye en sentido general, pero es tan débil la relación que no nos da ninguna información relevante.
Ahora vemos esta relación graficada:
Mostrar código
edad_data <-data.frame(edad =sort(unique(train$edad)),mean_numerosi =tapply(train$numerosi, train$edad, mean))ggplot(edad_data, aes(x = edad, y = mean_numerosi)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5) +theme_minimal() +labs(title ="Relación entre Edad y Número de Siniestros",x ="Edad",y ="Promedio de Número de Siniestros" ) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
En efecto, no queda muy clara la relación entre las variables. Por lo tanto, tal como hemos aprendido, discretizar la variable puede ayudar a su interpretación. La discretizamos usando la función que habíamos creado previamente.
edad_data_discretizada <-data.frame(edad_discretizada =levels(train$edad_discretizada),mean_numerosi =tapply(train$numerosi, train$edad_discretizada, mean))edad_data_discretizada$edad_discretizada <-factor(edad_data_discretizada$edad_discretizada, levels = edad_data_discretizada$edad_discretizada)ggplot(edad_data_discretizada, aes(x =as.numeric(edad_discretizada), y = mean_numerosi)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5, aes(group =1)) +scale_x_continuous(breaks =1:length(edad_data_discretizada$edad_discretizada), labels = edad_data_discretizada$edad_discretizada) +theme_minimal() +labs(title ="Relación entre Edad Discretizada y Promedio de Número de Siniestros",x ="Grupo de Edad Discretizado",y ="Promedio de Número de Siniestros" ) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
Ahora la relación queda mucho más clara. Los conductores más jóvenes (36-40 años) tienen un promedio de siniestros bastante alto, probablemente debido a la falta de experiencia al volante y comportamientos de conducción más arriesgados. Los conductores de edad joven-media (20-36) muestra una disminución en el promedio de número de siniestros, lo cual indica un mejor manejo al volante o hábitos más seguros al volante con los años. En los conductores de edad media (40-53) disminuye el promedio de siniestros, lo cual puede deberse a que los conductores en este grupo de edad, además de tener más experiencia al volante, a menudo tienen una vida profesional y familiar activa, por lo que podría significar que conducen con más cuidado. Los conductores más mayores (53-84) tienen un promedio de siniestros más bajo, posiblemente porque conducen con mayor precaución o debido a que no suelen conducir largos y complicados viajes sino cortos viajes del día a día para lo esencial.
3.2.7 Sexo conductor
En este caso analizamos la variable del sexo del conductor. Esta variable vuelve a ser de tipo factor así que la vamos a analizar de forma parecida a la variable de tipo ed vehículo.
table(train$sexocondu)
M V
131 536
De acuerdo con los datos, la mayoría de las conductoras son femeninas (dando por hecho que M es masculino y V femenino), ya que hay 536 observaciones de conductores femeninos en comparación con 131 de conductores masculinos. Esto indica que el grupo de conductores en este conjunto de datos está desbalanceado en favor de las mujeres.
tapply(train$numerosi, train$sexocondu, mean)
M V
0.5648855 0.4197761
Como vemos, los hombres en esta base de datos tienden a tener más siniestros en comparación con las mujeres, pero debemos comprobar si hay diferencias significativas en el número de siniestros entre estos. Usando ANOVA sabemos si el sexo tiene un impacto importante en el número de siniestros o si es más probable que las diferencias que vemos sean solo por casualidad.
Analysis of Variance Table
Response: train$numerosi
Df Sum Sq Mean Sq F value Pr(>F)
train$sexocondu 1 2.22 2.21667 2.3671 0.1244
Residuals 665 622.75 0.93646
Recordando los umbrales comunes de significancia, el análisis de varianza muestra que el sexo tiene un impacto significativo al nivel de significatividad 0.01 en el número de siniestros. El valor p es bajo (0.1244), lo que sugiere que las diferencias en el número de siniestros entre los dos sexos podrían deberse. La variable sexo del conductor parece ser un factor importante para predecir el número de siniestros.
Este intervalo de confianza sugiere que con un 95% de confianza, el efecto de ser mujer en la variable objetivo es negativo. Es decir, el ser mujer reduce el número de siniestros en comparación con los hombres tal como esperábamos con los estudios anteriores.
3.2.8 Antigüedad en la compañía
Ahora estudiamos la variable antigüedad en la compañía:
Mostrar código
ggplot(train, aes(x = anticia)) +geom_density(fill = colores[3], color = colores[5], alpha =0.6, bw =0.5) +theme_minimal() +labs(title ="Distribución de la Antigüedad en la Compañía",x ="Antigüedad del Vehículo (años)",y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
La mayoría de los clientes en la muestra llevan entre 0 y 5 años en la compañía Este rango tiene el pico de densidad más alto, cerca de 0.20, lo que sugiere que los clientes nuevos son los más comunes. A medida que la antigüedad aumenta, la densidad disminuye gradualmente. Esta distribución nos recuerda a una de tipo exponencial. Las vemos comparadas a continuación:
Como podemos observar, las distribuciones son bastante similares. Sin embargo, la distribución real no comienza tan alta como la exponencial y presenta un descenso más pronunciado entre aproximadamente los 3 y los 7 años. Además, hacia el final, la distribución real no muestra una caída tan pronunciada, sino más mantenida. Aunque las dos distribuciones son similares, las diferencias en el inicio, el comportamiento intermedio y el final indican que la distribución real tiene características particulares que no se ajustan perfectamente al modelo exponencial, sugiriendo que otros factores están influyendo en la antigüedad de los clientes en la compañía.
Ahora estudiamos su relación con la variable respuesta:
cor(train$anticia, train$numerosi)
[1] -0.01805373
Observamos una correlación muy débil (\(r=-0.01\)), lo que sugiere que, en general, a medida que un cliente lleva más años con la empresa, tiende a tener menos accidentes con el vehículo. Sin embargo, esta relación es tan débil que no es lo suficientemente “significativa” como para considerarse relevante.
Mostrar código
anticia_data <-data.frame(anticia =sort(unique(train$anticia)),mean_numerosi =tapply(train$numerosi, train$anticia, mean))ggplot(anticia_data, aes(x = anticia, y = mean_numerosi)) +geom_point(color = colores[5], size =2) +# Puntos en el gráficogeom_line(color = colores[3], size =0.7, aes(group =1)) +# Línea de conexión entre puntoslabs(title ="Relación entre Antigüedad del Vehículo y Promedio de Siniestros",x ="Antigüedad del Vehículo (años)",y ="Promedio de Número de Siniestros" ) +theme_minimal() +# Estilo minimalistatheme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )
De nuevo, el gráfico no deja clara ninguna relación entre las variables, así que discretizamos la variable agrupándola en si lleva menos de 3 años en la empresa (nuevo o new) o más (veterano o vet).
A pesar de que los vehículos nuevos tienen un promedio ligeramente mayor de siniestros que los vehículos veteranos, la diferencia es bastante pequeña. Esto sugiere que la antigüedad del vehículo, según esta discretización en “New” y “Vet”, no parece tener un impacto significativo en el número promedio de siniestros.
3.2.9 Usos del Vehículo
Por último, analizamos la variable uso del vehículo que ya habíamos discretizado en dos grupos al inicio de la práctica.
table(train$usovehi2)
G1 G2
627 40
tapply(train$numerosi, train$usovehi2, mean)
G1 G2
0.4529506 0.3750000
Los vehículos con un uso más alto (grupo G1 con uso mayor a 500) tienen un promedio de siniestros más alto en comparación con los vehículos con un uso menor (grupo G2 con uso menor a 100). Esto podría sugerir que los vehículos con más uso (G1) están involucrados en más accidentes o siniestros en promedio.
La convertimos a tipo factor para hacer un análisis más profundo de estas diferencias:
Analysis of Variance Table
Response: train$numerosi
Df Sum Sq Mean Sq F value Pr(>F)
train$usovehi2 1 0.23 0.22848 0.2432 0.6221
Residuals 665 624.74 0.93945
Dado que el valor p es alto, 0.6221, no hay evidencia estadística suficiente para concluir que el uso del vehículo tiene un efecto significativo sobre el número de siniestros. En otras palabras, las diferencias entre los grupos G1 y G2 en términos de su promedio de siniestros podrían darse simplemente por el azar.
El intervalo de confianza que incluye el valor 0 indica que la diferencia en el efecto entre el grupo G1 y el grupo G2 no es estadísticamente significativa. Por lo tanto, como ya hemos visto, no podemos afirmar que ser parte del grupo G2 tenga un impacto significativo en el número de siniestros.
3.3 Algunas conclusiones del EDA
En el análisis exploratorio de datos (EDA), la mayoría de las variables estudiadas y sus respectivas transformaciones parecen tener una capacidad predictiva notable por sí solas, con algunas excepciones como las variables plazas del vehículo y tipo de vehículo, que no muestran un impacto significativo. Por tanto, se optará por incluir todas estas variables en el modelo inicial.
En el caso de edad o antigüedad del carnet, hemos observado una relación no lineal con la variable dependiente. Este comportamiento sugiere que la edad podría no seguir una tendencia lineal simple, y es por eso que hemos decidido incluir transformaciones de la edad, como edad al cuadrado y edad al cubo, que pueden captar esta no linealidad de manera más efectiva, que funciona como alternativa al uso de los modelos GAM.
# Edad cuadrado y cubotrain$edad2 <- train$edad^2train$edad3 <- train$edad^3
Hay muchas oportunidades para construir más predictores que podrían mejorar el modelo. Algunas posibles ideas incluyen:
Discretizar variables adicionales o hacerlo de manera diferente, por ejemplo, agrupando rangos de edad, antigüedad del vehículo, etc.
Incorporar efectos de interacción entre variables, ya que algunas combinaciones de variables podrían tener un impacto mayor al combinarse.
Crear nuevas variables derivadas a partir de las que ya tenemos disponibles. Un ejemplo sería utilizar el código postal para generar una nueva variable que capture la localización geográfica y su influencia potencial en el comportamiento de los siniestros o cualquier otra variable de interés.
Esto puede enriquecer el modelo y mejorar su capacidad predictiva, por lo que seguir explorando nuevas combinaciones y transformaciones de variables sería una parte clave del proceso de modelado en otro caso.
4 Evaluando modelos
Propuestas de modelización
Especificación:
Variables que van a entrar en el predictor lineal
Distribución para la variable respuesta
Función link
Y con enfoque ML, tipo de regularización/penalización
Distribuciones para la variable respuesta
Como modelos para la variable respuesta podemos testar:
Poisson
Binomial Negativa
ZI (inflada en cero) + Poisson
ZI (inflada en cero) + BN
Normal
Funciones link
Como funciones link:
Logaritmo
Identidad
Regularización
Penalización: * Net Elastic (Lasso, Ridge)
Modelos sin y con regularización
Modelos sin regularización: Estos modelos ajustan los datos de forma directa, maximizando el ajuste sin restricciones. Aunque pueden ofrecer un buen desempeño inicial, son más propensos al sobreajuste, capturando ruido en lugar de patrones generales.
Modelos con regularización: Incorporan penalizaciones que limitan la complejidad del modelo, favoreciendo soluciones más simples y generalizables. Ayudan a evitar el sobreajuste y suelen ser más robustos en escenarios reales..
4.1 Modelos sin regularización
Cuando utilizamos modelos sin regularización hemos de estar atentos a los problemas de multicolinealidad.
Antes de comenzar a modelizar hemos de eliminar algunas variables por que es lo razonable, lógico, desde el punto de vista de modelización y de computación.
Eliminamos codigopostal y ccaa ya que no aportan información relevante al modelo y no necesariamente tienen una relación significativa o directa con el número de siniestros, ya que ambas variable son identificadores geográficos y podrían afectar a la complejidad del modelo en lugar de mejorarlo. También eliminamos las variables iniciop y finalp, debido a que ya hemos calculado con anterioridad la diferencia entre estas para determinar la duración de la póliza, por lo que mantener estas variables sería redundante y podría dificultar la interpretación del modelo.
Po.log <-glm(numerosi ~ . , data = train, family =poisson(link ="log"))# Obtenemos los valores iniciales y nos aseguramos que no haya coeficientes negativos o NAstarting.values <-coef(Po.log)starting.values[starting.values <0] <-0.0001starting.values[is.na(starting.values)] <-0.0001# Ajuste modelo Poisson con enlace identidad usando los valores inicialesPo.id <-glm(numerosi ~ . , data = train, start = starting.values,family =poisson(link ="identity"))
Ajuste modelo (clásico) Binomial Negativo
BN.log <-glm.nb(numerosi ~ . , data = train, link = log)# starting.values <- coef(BN.log)# starting.values[starting.values < 0] <- 0.0001# BN.id <- glm.nb(numerosi ~ . , data = train,# start = starting.values, link = identity)
Ajustes de modelos (clásicos) inflados de ceros
# Modelo inflado de ceros con PoissonZI.Po <-zeroinfl(numerosi ~ ., data = train, dist ="poisson")# Modelo inflado de ceros con Binomial NegativaZI.BN <-zeroinfl(numerosi ~ ., data = train, dist ="negbin")
Normal
La variable de respuesta, numerosi, representa el número de siniestros para cada observación, por lo que se trata de una variable discreta que toma valores enteros no negativos entre los cuales es común observar una gran cantidad de ceros.
Debido a esto, realizar un ajuste con modelo Normal no tiene sentido ya que la Normal es una distribución que admite valores negativos y podría generar errores significativos en la estimación de los parámetros y predicciones poco precisas. Por lo tanto, descartamos este modelo.
4.1.1 El rol del tiempo de exposición: Modelos con Offsets (y pesos)
A continuación, realizamos diversos modelos para analizar el número de siniestros, incluyendo el manejo del tiempo de exposición utilizando offsets y pesos.
Mostrar código
# Poisson con offsetPo.log.off <-glm(numerosi ~ . ,data =subset(train, select =-c(duracion_poliza)),offset =log(train$duracion_poliza /365.25), # Convertimos la duración a años family =poisson(link ="log"))# Poisson con pesos (corrigiendo `weight` a `weights`)Po.log.wgt <-glm(numerosi ~ . ,data =subset(train, select =-c(duracion_poliza)),weights = train$duracion_poliza /365.25, family =poisson(link ="log"))
El término offset nos permite incorporar la duración de la póliza como una variable que ajusta la tasa de siniestros por unidad de tiempo, garantizando que la duración de la póliza se considere de manera proporcional al análisis. En este caso, la duración se convierte a años, y su efecto se incluye de manera multiplicativa en el modelo de Poisson, ajustando las tasas de siniestros en función del tiempo de exposición al riesgo. Por otro lado, el modelo con pesos utiliza la duración de la póliza como un peso, asignando mayor relevancia a las observaciones con una mayor duración, y corrigiendo así el modelo para que tenga en cuenta de manera adecuada la variable de tiempo sin necesidad de usar un término offset explícitamente.
Observamos que el modelo con offset presenta un AIC más bajo que el modelo con pesos. Esto sugiere que el primer modelo se ajusta mejor a los datos, además de ser más interpretable en términos de tasas de siniestros ajustadas por la duración de la póliza.
Realizamos otros modelos utilizando una fórmula más detallada, especificando los predictores. La variable respuesta aparece como un offset explícito, lo cual aclara su contribución al modelo.
La Binomial Negativa permite manejar la sobredispersión que puede llegar a ocurrir con Poisson, y los modelos inflados de cero se utilizan cuando hay una gran cantidad de ceros en la variable dependiente, en este caso el número de siniestros, condición que no puede ser capturada adecuadamente por Poisson o Binomial Negativa.
4.1.2 Evaluación clásica de los modelos
Procedemos a comparar los diferentes modelos predictivos utilizando el criterio AIC.
AICs <-c(AIC(Po.log), AIC(Po.id), AIC(BN.log), AIC(ZI.Po), AIC(ZI.BN),AIC(Po.log.off), # ‘Ajuste’ por número de observacionesAIC(Po.log.wgt)*nrow(train)/sum(train$duracion_poliza/365.25),AIC(BN.log.off), AIC(ZI.Po.off), AIC(ZI.BN.off))
colores1 <-carto_pal(11, "Teal")barplot( errores_clasicos$AICs, names.arg = errores_clasicos$modelo, col = colores1[rank(errores_clasicos$AICs)], border ="white", main ="Valores AIC por Modelo", xlab ="Modelos", ylab ="AIC", las =1, cex.names =0.6,cex.axis =0.8, cex.main =1.2)
Vemos que el modelo Poisson inflado de ceros con offset presenta el AIC más bajo (AIC = 1098.081), lo cual nos indica que es el modelo que mejor se ajusta a los datos. Elegiríamos este modelo, ya que un AIC bajo sugiere que este tiene una mejor relación entre la capacidad de ajuste y el número de parámetros, evitando así el sobreajuste.
4.1.3 Evaluación enfoque ML
Definición de nuevos predictores definidos en el conjunto de training en test.
Discretizamos la potencia usando cuantiles en el test, como ya hicimos anteriormente en el training set. También clasificamos la antigüedad de la compañía en “Nueva” o “Veterana”, creamos tranformaciones para capturar relaciones no lineales entre edad y otras variables en modelos predictivos, eliminamos columnas innecesarias y por último, filtramos los datos para eliminar observaciones donde la duración de la póliza sea 0 para evitar errores posteriores.
Eliminación de observaciones
Tras haber intentado generar las predicciones, obtuvimos un error sobre los niveles 5 y 19 de la variable prov. En esta tabla observamos que dichos niveles no están presentes en el conjunto de entrenamiento mientras que sí lo están en el conjunto de prueba, generando un error ya que los modelos no sabrán tratar estos niveles. Para evitar errores durante la generación de predicciones y demás, eliminamos las observaciones de dichos niveles en el test set.
Pasamos a generar las predicciones en el conjunto de prueba para los modelos sin regularización. Descartamos el modelo Poisson con pesos debido a que, como hemos podido observar, tiene un AIC significativamente más alto respecto al resto de modelos.
El gráfico muestra la comparación de las distribuciones de las predicciones obtenidas con distintos modelos en el test set, frente a los datos reales del mismo. Observamos que en la mayoría de modelos las predicciones siguen una tendencia parecida a la de los datos reales, con un pico inicial elevado cerca de 0 sinestros. Mientras que modelos como Po_id presentan diferencias marcadas respecto a realidad, otros como ZI_BN_off parecen ajustarse mejor.
No obstante, para determinar correctamente cuál es el mejor modelo, procedemos a calcular los errores cuadráticos y absolutos.
El error cuadrático se utiliza para evaluar el ajuste de los modelos comparando las predicciones con los valores reales observados en el test set. Penaliza más fuertemente las grandes desviaciones entre las predicciones y los valores reales.
A diferencia del error cuadrático, el error absoluto no penaliza en exceso las grandes desviaciones, sino que simplemente calcula la suma de las diferencias absolutas entre las predicciones y los valores reales. Esto es útil cuando se quiere tener una idea más directa de la magnitud de los errores sin que los grandes errores afecten de forma desproporcionada el cálculo total.
Resumen errores (modelos sin regularización)
A continuación, se muestra una tabla resumiendo cuadráticos y absolutos obtenidos para cada modelo para facilitar la comparación entre los modelos ajustados.
Atendiendo a los errores y a las características de los modelos, elegiríamos el modelo de la Binomial Negativa. Como se ha comentado anteriormente, la distribución Binomial Negativa se ajusta mejor a nuestra muestra y es útil a la hora de sobrellevar la sobredispersión en los datos que se genera comúnmente en variables donde la varianza supera la media. A pesar de que los errores más bajos los presenta el modelo Poisson con enlace log, los valores tanto del error cuadrático como del absoluto son muy cercanos entre ambos modelos.
Respecto a la evaluación clásica, aunque evaluar el AIC puede resultar en modelos más complejos con menos sobreajuste, los errores nos ofrecen una visión directa del rendimiento predictivo del modelo. Esto puede hacer que un modelo con un AIC más bajo pero con peores errores no será necesariamente la mejor opción si nos fijamos en la precisión de las predicciones.
4.2 Modelos con regularización
Tras una comparación inicial de modelos sin regularización basándonos en el ajuste de los datos, ahora incluiremos la regularización, penalizando la complejidad del modelo y promoviendo soluciones más simples y robustas.
Los Paquetes glmnet y mpath
Cuestiones sobre las funciones - Conversión de data.frame en matriz (construyendo variables dummy para factores e intercepto), necesario para predecir.
model.matrix(~., data)
Modelo con respuesta Poisson
cv.glmnet(x, y, family=“poisson”, alpha, nfolds) glmnet(x, y, family=“poisson”, alpha) cvglmreg(formula, data, family=“poisson”, penalty, nfolds) glmreg(formula, data, family=“poisson”, penalty)
# Convertimos los data.frame en objetos matriz respetando factorestrain.m <-model.matrix(~., data=train)train.m <- train.m[,colnames(train.m) !="numerosi"]test.m <-model.matrix(~., data=test0)test.m <- test.m[,colnames(test.m) !="numerosi"]
Modelos con regularización
Se estima un modelo, existe una solución, por cada valor del parámetro de penalización.
Descarte de modelos
A lo largo del resto del trabajo, hemos decidido descartar algunos modelos debido a errores técnicos que no se pudieron solucionar durante el proceso.
4.2.1 Ajustes con glmnet
Mostrar código
# Respuesta NormalNo.re <-glmnet(x = train.m, y = train$numerosi, alpha=0.5)# Respuesta PoissonPo.re <-glmnet(x = train.m, y = train$numerosi,family ="poisson", alpha =0.5)# Respuesta Poisson con offsetPo.re.off <-glmnet(x = train.m, y = train$numerosi, alpha =0.5, offset =log(train$duracion_poliza/365.25), family ="poisson")
Utilizamos glmnet para ajustar modelos con regularización mediante la penalización Elastic Net (alpha = 0.5), que combina las propiedades de las penalizaciones Lasso y Ridge.
4.2.2 Ajustes con mpath
Mostrar código
# Respuesta PoissonPo.re2 <-glmreg(numerosi ~ . , data = train,family ="poisson", penalty="enet")# Respuesta Binomial NegativaBN.re <-glmregNB(numerosi ~ . , data = train, penalty="enet")# Respuestas Infladas de CerosZI.Po.re <-zipath(numerosi ~ .| . , data= train,family ="poisson", link ="logit", penalty="enet")ZI.BN.re <-zipath(numerosi ~ .| . , data= train,family ="negbin", link ="logit", penalty="enet")# Modelos Po y BN con offsetformula3 <- numerosi ~ potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia + g.anticia + edad2 + edad3# Modelos con offsetPo.re2.off <-glmreg(formula3, data = train, family ="poisson",offset =log(duracion_poliza/365.25), penalty="enet")BN.re.off <-glmregNB(formula3, data = train, offset =log(duracion_poliza/365.25), penalty="enet")# Modelos ZI con offsetformula4 <-formula(numerosi ~ potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia + g.anticia + edad2 + edad3 | potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia + g.anticia + edad2 + edad3)ZI.Po.re.off <-zipath(formula4, data= train, family ="poisson", offset =log(duracion_poliza/365.25), link ="logit", penalty="enet", maxit =100)ZI.BN.re.off <-zipath(formula4, data= train, family ="negbin", offset =log(duracion_poliza/365.25), link ="logit", penalty="enet", maxit =100)
Ajustamos los modelos con mpath con fórmulas específicas para asegurarnos de que las variables predictoras relevantes estén correctamente incluidas en los ajustes. Este enfoque nos permite comparar diferentes técnicas y distribuciones, considerando características particulares de los datos y explorando la influencia de la regularización en cada caso.
4.2.3 Elección del modelo
Disponemos de un modelo para cada lambda testado, ¿cuál elegimos?
Vamos a considerar dos alternativas:
Predecimos con el modelo con el lambda que genera el “mejor” ajuste con medida clásica (AIC).
Elegimos el mejor modelo por Cross-Validation
Elección de modelos con regularización basándonos en el AIC y predicción
Mostrar código
# Modelo de respuesta Normalselec <-which(No.re$dev.ratio ==max(No.re$dev.ratio))p.No.re <-predict(No.re, test.m, type ="response")[, selec]# Modelo de respuesta Poisson (glmnet)selec <-which.max(Po.re$dev.ratio)p.Po.re <-predict(Po.re, test.m, type ="response")[, selec]# Modelo de respuesta Poisson y offset glmnetselec <-which.max(Po.re.off$dev.ratio)p.Po.re.off <-predict(Po.re.off, test.m, newoffset =log(test$duracion_poliza /365.25),type ="response")[, selec]# Modelo de respuesta Poisson (mpath)selec <-which.min(Po.re2$resdev)p.Po.re2 <-predict(Po.re2, test0, which = selec, type ="response")# Modelo de respuesta Poisson y offset (mpath)selec <-which.min(Po.re2.off$resdev)p.Po.re2.off <-exp(predict(Po.re2.off, test0, which = selec,newoffset =log(test0$duracion_poliza /365.25)))# Modelo de respuesta BN (mpath)selec <-which.min(AIC(BN.re))p.BN.re <-predict(BN.re, test0, which = selec, type ="response")# Modelo de respuesta BN y offset (mpath)selec <-which.min(AIC(BN.re.off))p.BN.re.off <-exp(predict(BN.re.off, test0, which = selec,newoffset =log(test0$duracion_poliza /365.25)))# Modelo de respuesta ZI-Po (mpath)selec <-which.min(AIC(ZI.Po.re))p.ZI.Po.re <-predict(ZI.Po.re, test0, which = selec,type ="response")# Modelo de respuesta ZI-Po y offset (mpath)selec <-which.min(AIC(ZI.Po.re.off))p.ZI.Po.re.off <-predict(ZI.Po.re.off, test0, which = selec,type ="response")# Modelo de respuesta ZI-BN (mpath)selec <-which.min(AIC(ZI.BN.re))p.ZI.BN.re <-predict(ZI.BN.re, test0, which = selec, type ="response")# Modelo de respuesta ZI-BN y offset (mpath)selec <-which.min(AIC(ZI.BN.re.off))p.ZI.BN.re.off <-predict(ZI.BN.re.off, test0, which = selec, type ="response")
Realizamos la selección y predicción de varios modelos ajustados con regularización, considerando qué tanto del comportamiento de los datos explica el modelo, o qué tan bueno es el modelo según el criterio del AIC. Además, incorporamos ajustes con offset para manejar correctamente la exposición al riesgo, asegurando que los modelos sean más interpretables y relevantes para el contexto.
Comparamos las densidades de los predicciones de varios modelos con las de los datos reales del número de siniestros. Los modelos BN_re y BN_re_off destacan por ajustarse mejor a los datos reales, especialmente en los valores más frecuentes (0 y 1 siniestros). En cambio, modelos como Po_re_off tienden a sobreestimar los siniestros altos, lo que indica un ajuste menos preciso. Por otro lado, el modelo No_re, aunque tiene bajos errores, subestima frecuencias importantes y no representa bien la variabilidad de los datos. En general, los modelos basados en la Binomial Negativa logran capturar mejor las características reales de la distribución.
Ahora, al igual que hicimos antes con los modelos sin regularización, comprobamos los errores cuadráticos y absolutos.
errores_combinados <-rbind( errores_re_AIC$errores_cuadraticos_re_AIC, errores_re_AIC$errores_abs_re_AIC)colnames(errores_combinados) <-c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", "BN.re", "ZI.Po.re", "Po.re2.off","BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off")barplot(errores_combinados, beside =TRUE, names.arg =c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", "BN.re", "ZI.Po.re", "Po.re2.off","BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off"), col =c(colores2[2], colores2[7]),border ="white", main ="Comparación de Errores Cuadráticos y Absolutos", xlab ="Modelo", ylab ="Error", las =1, cex.names =0.6,cex.axis =0.8, cex.main =1.2)
Elegiríamos el modelo de la Binomial Negativa con offset. Entre los modelos evaluados, este modelo es el que presenta los mejores resultados generales, con el mínimo error tanto absoluto (163.79), como cuadrático (213.78). El valor del error cuadrático indica que el modelo maneja eficazmente posibles desviaciones extremas en los datos, mientras que el error absoluto refleja que en promedio sus predicciones se acercan más a los valores reales comparado con el resto de modelos.
Comparando con los modelos sin regularización, observamos que los errores de los modelos con regularización son más bajos en general. Existe una disminución de los errores en el modelo regularizado que hemos elegido respecto al modelo sin regularización anteriormente preferido.
4.2.4 Elección de modelos usando cross-validation
Mostrar código
train.m <-model.matrix(~., data=train)train.m <- train.m[,colnames(train.m) !="numerosi"]test.m <-model.matrix(~., data=test0)test.m <- test.m[,colnames(test.m) !="numerosi"]# Convertimos a matriz dispersatrain.m <-as(train.m, "dgCMatrix")# Modelo de respuesta Normalcv.No <-cv.glmnet(x=train.m, y=train$numerosi, alpha=0.5)# Modelo de respuesta Poisson (glmnet)cv.Po <-cv.glmnet(x=train.m, y=train$numerosi, alpha=0.5, family ="poisson")# Modelo de respuesta Poisson y offset (glmnet)cv.Po.off <-cv.glmnet(x = train.m, y = train$numerosi, alpha =0.5, family ="poisson", offset =log(train$duracion_poliza/365.25))# Modelo de respuesta Poisson (mpath)#cv.Po2 <- cv.glmreg(numerosi ~ . , data = train, family="poisson",# penalty = "enet")# Modelo de respuesta Poisson y offset (mpath)# Ajustar el modelo Poisson con penalización#cv.Po2.off <- cv.glmreg(formula3 , data = train_log, # family = "poisson", penalty="enet",# offset = log(duracion_poliza/365.25))# Modelo de respuesta BN (mpath)#cv.BN <- cv.glmregNB(formula3 , data= train, # penalty="enet", nfolds = 3, trace = TRUE)# Modelo de respuesta BN y offset (mpath)#cv.BN.off <- cv.glmregNB(formula3, data = train, # offset = log(train$duracion_poliza/365.25))# Modelo de respuesta ZI-Po (mpath)#cv.ZI.Po <- cv.zipath(numerosi ~ .| . , data= train,# family="poisson", link = "logit",# penalty="enet")# Modelo de respuesta ZI-BN (mpath)#cv.ZI.BN <- cv.zipath(numerosi ~ .| . , data= train, # family="negbin", link = "logit",# penalty="enet")# Modelo de respuesta ZI-Po y offset (mpath)# Variables y objetos auxiliarestrain$duracion_poliza.a <- train$duracion_poliza/365.25test0$duracion_poliza.a <- test0$duracion_poliza/365.25formula5 <-formula(numerosi~ potencia + antivehi+ tipovehi+ plazasvehi+ edad + sexocondu+ anticia+ prov+ usovehi2 + g.potencia+ g.anticia+ edad2 + edad3 + edad_discretizada+offset(log(duracion_poliza.a))| potencia + antivehi+ tipovehi+ plazasvehi+ edad + sexocondu+ anticia+ prov+ usovehi2 + g.potencia+ g.anticia+ edad2 + edad3 + edad_discretizada+offset(log(duracion_poliza.a)))#cv.ZI.Po.off<-cv.zipath(formula5, data= train, link = "logit",# family= "poisson", penalty="enet")# Modelo de respuesta ZI-BN y offset (mpath)#cv.BN.Po.off <- cv.zipath(formula5, data= train, family= "negbin",# link = "logit", penalty="enet")
El enfoque de cross-validation nos permite comparar el rendimiento de los modelos en diferentes conjuntos de entrenamiento y prueba, ayudando a identificar el modelo más adecuado según las métricas de error generadas durante la validación.
Tras entrenar los modelos utilizando técnicas de selección de hiperparámetros mediante validación cruzada, generamos las predicciones de los mismos.
Mostrar código
resultados <-data.frame(reales = test0$numerosi,cv_No = p.cv.No,cv_Po = p.cv.Po,cv_Po_off = p.cv.Po.off,cv_ZI_Po_off = p.cv.ZI.Po.off,cv_ZI_BN_off = p.cv.ZI.BN.off)resultados <-melt(resultados, id.vars ="reales", variable.name ="Modelo", value.name ="Predicciones")resultados$Modelo <-recode(resultados$Modelo,"lambda.1se"="cv_No","lambda.1se.1"="cv_Po","lambda.1se.2"="cv_Po_off")ggplot(resultados, aes(x = Predicciones, fill = Modelo)) +geom_density(alpha =0.4) +geom_density(data = resultados, aes(x = reales), fill ="black", alpha =0.2) +facet_wrap(~ Modelo) +scale_x_continuous(limits =c(0, 5)) +scale_fill_manual(values = colores3) +labs(x ="Número de Siniestros", y ="Densidad", title ="Comparación de Predicciones con Datos Reales") +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),axis.text.x =element_text(),panel.border =element_blank(),plot.margin =margin(2, 1, 1, 1, "cm"),legend.position ="none")
Al igual que antes, el gráfico compara la densidad de las predicciones, esta vez de los modelos validados con validación cruzada. El modelo Binomial con offset muestra el mejor ajuste, logrando capturar de forma precisa la distribución de los datos reales, especialmente en los siniestros más frecuentes. En contraste, cv_No subestima los valores intermedios, y cv_Po_off tiende a sobreestimar siniestros altos, reduciendo su precisión. Los resultados refuerzan la preferencia por cv_ZI_BN_off, que demuestra una mejor capacidad de generalización.
Volvemos a elegir el modelo BN.re.off como la mejor opción. Aunque fijándonos únicamente en el valor de los errores el modelo No.re podría ser más adecuado, este no incluye mecanismos para manejar sobredispersión ni incorpora términos de regularización que ayuden a evitar el sobreajuste
El modelo de la Binomial Negativa es particularmente adecuado para nuestro conjunto de datos, ya que, como se ha mencionado, maneja eficientemente la sobredispersión que caracteriza a la variable dependiente. Además, la regularización, a través del término offset, contribuye a evitar problemas de sobreajuste y mejora la estabilidad del modelo.
Comparando con los modelos anteriores, BN.re.off demuestra una mejora significativa respecto a los modelos sin regularización, como BN.log.off, que aunque presenta buenos resultados, no alcanza el desempeño del modelo seleccionado en las métricas clave. En relación a los modelos evaluados únicamente con AIC (evaluación clásica), BN.re.off logra un balance superior entre ajuste y capacidad predictiva, superando a opciones como ZI.Po.off y ZI.BN, cuyos valores de AIC más bajos no se traducen en un mejor desempeño en validación cruzada.
5. Conclusiones
En conclusión, dado que nuestro estudio se basa en una muestra donde la varianza de los datos supera la media, la Binomial Negativa es preferible gracias a su capacidad para capturar la estructura de estos datos, manejar la sobredispersión y garantizar predicciones confiables en nuevos conjuntos.
La decisión de reducir la muestra nos permitió agilizar la convergencia de los modelos y continuar con el análisis. Sin embargo, es importante destacar que esto podría haber limitado la capacidad de los modelos para captar patrones más complejos presentes en la muestra completa. Adicionalmente, se eliminaron algunas observaciones específicas para evitar problemas al realizar las predicciones, como aquellas correspondientes a niveles de factores no presentes en el conjunto de entrenamiento. Estos ajustes fueron necesarios para garantizar la ejecución, pero podrían haber afectado la representatividad de los resultados. Por lo tanto, sería interesante repetir el análisis utilizando la totalidad de los datos originales en un equipo de soporte más preparado para manejar tal cantidad de datos.
[INE - Instituto Nacional de Estadística. (n.d.). INEbase/ Clasificaciones /Relación de municipios, provincias, comunidades y ciudades autónomas con sus códigos / Relación de comunidades y ciudades autónomas con sus códigos. INE - Instituto Nacional De Estadística.] (https://www.ine.es/daco/daco42/codmun/cod_ccaa_provincia.htm)
[¿Qué es la correlación estadística y cómo interpretarla? (n.d.). Máxima Formación.] (https://www.maximaformacion.es/blog-dat/que-es-la-correlacion-estadistica-y-como-interpretarla/)
[colaboradores de Wikipedia. (2023, November 22). Distribución log-normal. Wikipedia, La Enciclopedia Libre.] (https://es.wikipedia.org/wiki/Distribuci%C3%B3n_log-normal)
[Color palettes – Data Visualization with R. (n.d.). Data Visualization With R.] (https://datavizs24.classes.andrewheiss.com/resource/colors.html)
[What are continuous probability distributions & their 8 common types? | KNIME. (n.d.). KNIME.] (https://www.knime.com/blog/continuous-probability-distribution)
Source Code
---title: "Pricing II: Aplicación de R a la tarificación"authors: "Carolina Armella, Ágatha del Olmo y María Yu García"date: "2024-10-18"format: html: embed-resources: true theme: zephyr code-tools: true toc: true toc-depth: 4 execute: warning: false---```{r setup, include=FALSE}knitr::opts_chunk$set( warning = FALSE, message = FALSE, comment = "", fig.align = "center", fig.show = "hold", fig.height = 4, fig.width = 8, out.width = "80%")``````{r librerias, warning=FALSE, include=FALSE}library(rcartocolor)library(kableExtra)library(tidyverse)library(vcd)library(MASS)library(pscl)library(glmnet)library(mpath)library(corrplot)library(dplyr)library(ggplot2)library(reshape2)library(caret)colores <- carto_pal(7, "Teal")colores2 <- carto_pal(7, "RedOr")```## 1. IntroducciónEl propósito de este análisis es predecir el número esperado de reclamaciones que realizará una póliza a partir de la experiencia de siniestralidad histórica de una cartera de seguros. Este estudio se basa en la implementación de modelos lineales generalizados (GLM), que permiten realizar predicciones precisas manteniendo un alto grado de interpretabilidad.Para este caso, utilizamos un conjunto de datos reales proporcionado por una compañía de seguros, que recoge información sobre 30,000 pólizas. En el presente análisis, para simplificar y reducir el tiempo de procesamiento, trabajamos con una muestra aleatoria de la base de datos original.El análisis incluye una exploración detallada de los datos (Exploratory Data Analysis o EDA) y la implementación de modelos de aprendizaje automático para evaluar diferentes metodologías. Además, dado el contexto regulatorio de la industria de seguros, donde la transparencia y la explicabilidad son fundamentales, se priorizan los modelos que cumplen con estos requisitos, por eso no vamos a usar modelos de tipo caja negra.## 2. MetodologíaPara implementar un sistema automatizado de evaluación crediticia (credit scoring), es fundamental contar con modelos que permitan identificar con precisión los riesgos asociados a los préstamos. Además, la normativa bancaria en muchos países, incluido España, exige que las instituciones expliquen a los solicitantes los motivos para aceptar o rechazar una solicitud de crédito. Por ello, los modelos empleados deben ser transparentes. Herramientas más avanzadas como las redes neuronales o los modelos de máquinas de soporte vectorial (SVM) son consideradas "cajas negras" y no son fáciles (y a veces son imposibles) de interpretar y por lo tanto no cumplen con este requisito de transparencia.### 2.1 Exploración de los datosLa base de datos corresponde a una extracción aleatoria de la BBDD de experiencia de siniestralidad de una compañía de seguros y consta de 19 variables capturando distintas características de las pólizas.::: {.callout-warning title="Reducción del dataset"}Hemos decidido reducir la base de datos a 1000 observaciones en vez de las 27378 originales para agilizar la convergencia por problemas de rendimiento.:::```{r read}datos <- read.csv("D:\\3 curso\\TP\\Grupo 2.csv", sep=";")datos <- datos[1:1000,] # reducimos el training set para que tarde menosdim(datos)str(datos)```En este caso, nuestra variable respuesta es numerosi, que indica el número de siniestros reportados por cada póliza. Esta variable es crucial porque el objetivo del modelo es predecir el número esperado de reclamaciones en función de las características de la póliza y del asegurado.La variable nploiza es simplemente un identificador único de cada póliza, y no contiene información relevante o explicativa para la predicción de siniestros, por lo tanto la podemos eliminar.Además, hay dos variables que nos dan información adicional. Podemos sustituir iniciop y finalp por la duración de la póliza en días o meses. Este predictor podría ser relevante teniendo en cuenta la exposición al riesgo. En cuanto al código postal, aunque en apariencia es solo una etiqueta numérica, contiene información geográfica valiosa, y para hacer más fácilmente interpretable este predictor extraemos la comunidad autónoma ya que no contamos con los metadatos de ccaa para saber a qué referencia cada número, por lo tanto esa variable la eliminamos.### 2.2 Preprocesamiento de los datosEliminamos por las razones indicadas la variable npoliza.```{r var}# Descartamos el número de pólizadatos <- subset(datos, select = -c(npoliza))``````{r}#| code-fold: true#| code-summary: "Mostrar código"datos <- datos %>%mutate(region =substr(as.character(codigopostal), 1, 2)) %>%mutate(comunidad =case_when( region %in%c("01", "48", "20") ~"País Vasco", region %in%c("31") ~"Navarra", region %in%c("03", "12", "46") ~"Comunidad Valenciana", region %in%c("21", "41", "14", "29", "11", "18", "23", "04") ~"Andalucía", region %in%c("08", "17", "25", "43") ~"Cataluña", region %in%c("15", "27", "32", "36") ~"Galicia", region %in%c("33") ~"Asturias", region %in%c("39") ~"Cantabria", region %in%c("09", "42", "47", "34", "37", "05", "24", "40", "49") ~"Castilla y León", region %in%c("45", "13", "16", "19", "02") ~"Castilla-La Mancha", region %in%c("22", "44", "50") ~"Aragón", region %in%c("30") ~"Murcia", region %in%c("28") ~"Madrid", region %in%c("07") ~"Islas Baleares", region %in%c("35", "38") ~"Canarias", region %in%c("06", "10") ~"Extremadura", region %in%c("26") ~"La Rioja", region %in%c("51", "52") ~"Ceuta y Melilla" ))datos <- datos[, -19]```Comprobamos que no haya NULL en la variable creada ya que hemos contemplado todas las posibles combinaciones de códigos postales y sus comunidades autónomas según el INE en la fórmula.```{r}sum(is.na(datos$comunidad))```Como vemos al crear la nueva variable comunidad autónoma, nos encontramos con muchos NA, por lo tanto, o la base de datos es completamente inventada o no sigue el sistema de códigos postales oficial, lo cual hace imposible su uso. Por lo tanto, eliminamos la variable creada.```{r}datos <- datos[, -19]```Muchas de las variables no tienen la clase que les correspondería. Primero las variables iniciop y finalp las pasamos a formato fecha. También creamos la variable duración póliza con la diferencia entre estas dos, que nos indica cuánto ha durado la póliza en días.```{r tables}# Convertimos en formato fecha: as.Date(x, "%m/%d/%Y") datos$iniciop <- paste0(substr(datos$iniciop, 5, 6), "/", substr(datos$iniciop, 7, 8), "/", substr(datos$iniciop, 1, 4) )datos$iniciop <- as.Date(datos$iniciop , "%m/%d/%Y")datos$finalp <- paste0(substr(datos$ finalp, 5, 6), "/", substr(datos$ finalp, 7, 8), "/", substr(datos$ finalp, 1, 4) )datos$finalp <- as.Date(datos$ finalp , "%m/%d/%Y")datos$duracion_poliza <- as.numeric(datos$finalp - datos$iniciop)```A veces también es necesario convertir variables a tipo factor. Los factores se usan comúnmente para representar variables categóricas, es decir, variables que toman un número limitado de valores. En nuestro caso:- tipovehi: Representar diferentes tipos de vehículos, en concreto 7.- codigopostal: Se convierte en un factor, ya que representa códigos postales, que son categorías y no valores numéricos continuos.- prov: Similar, representa provincias que se pueden tratar como categorías.- ccaa: Representa la comunidad autónoma (en España), lo que también es una variable categórica.```{r}datos$tipovehi <-as.factor(datos$tipovehi)datos$codigopostal <-as.factor(datos$codigopostal)datos$prov <-as.factor(datos$prov)datos$ccaa <-as.factor(datos$ccaa)str(datos)```Echamos un vistazo a las tablas de las variables no continuas para ver su frecuencia:```{r, include=FALSE}# Calculamos algunos resumenesapply(datos[,c(3, 5:11, 13, 15:18)], 2, table)```Creamos ahora una nueva variable llamada usovehi2 que clasifica la variable numérica usovehi en dos grupos, los que tienen un valor menor a 2000 que se pasan a llamar G1 y los que tienen un valor mayor, que se llamarán G2. Después eliminamos la variable original.```{r}datos$usovehi2 <-ifelse(datos$usovehi <20000, "G1", "G2")# Descartamos uso de vehiculo, poca varianzadatos <-subset(datos, select =-c(usovehi))```Observamos la matriz de correlaciones, que nos indica el coeficiente de correlación entre todos los pares de variables numéricas de la base de datos.::: {.callout-note title="Coeficiente de correlación"}El coeficiente de correlación mide la intensidad y dirección de la relación lineal entre dos variables:- $r=1$: Correlación positiva perfecta, a medida que una variable aumenta, la otra también aumenta proporcionalmente.- $r=−1$: Correlación negativa perfecta, a medida que una variable aumenta, la otra disminuye proporcionalmente.- $r=0$: No hay correlación lineal, las variables no tienen una relación lineal directamente.:::```{r}numericas <- datos[sapply(datos, is.numeric)]# Calcular la matriz de correlacióncorrelaciones <-cor(numericas, use ="pairwise.complete.obs")# Crear un gráfico de correlacióncorrplot(correlaciones, method ="circle", type ="upper", tl.col ="black", tl.srt =45, col = colores)```Como podemos ver, lógicamente las variables están perfectamente correlacionadas de forma directa con ellas mismas. Pero podemos ver otras relaciones interesantes, para empezar hay una correlación bastante alta entre la variable objetivo numerosi y las variables cuantiasi y siniestro. Descartamos siniestro y cuantiasi como predictores.```{r}cor(datos$siniestro, datos$numerosi)cor(datos$cuantiasi, datos$numerosi)datos <-subset(datos, select =-c(cuantiasi, siniestro))```Esto se debe a que la variable siniestro es un resultado de un evento ya ocurrido, y utilizarla para predecir el número de siniestros futuros es como usar la variable respuesta en el modelo de entrenamiento, lo que resulta en un modelo irrealista o sobreajustado por redundancias. Esto lo acabamos de ver reflejado en la matriz de correlaciones ($r=0.76$).Además, al estar cuantiasi también altamente correlacionada con la variable de respuesta ($r=0.35$), usarla como predictor puede causar multicolinealidad por redundancia, lo que afecta negativamente la interpretación y precisión de los coeficientes del modelo.El gráfico de correlaciones también nos ha mostrado una muy alta correlación entre potencia y valorvehi:```{r}cor(datos$valorvehi, datos$potencia)```Esto nos lleva a un alto riesgo de multicolinealidad, que ocurre cuando dos o más variables predictoras están altamente correlacionadas entre sí, lo que puede causar problemas en los modelos estadísticos, especialmente en los modelos de regresión, ya que puede dificultar la interpretación de los coeficientes tal como hemos mencionado con cuantiasi. Por lo tanto, eliminamos la variable de menos interés interpretativo para el estudio, valorvehi, ya que la potencia es algo más objetivo que no depende de la economía.```{r}datos <-subset(datos, select =-valorvehi)```Por último, el gráfico ha indicado una alta correlación entre edad y anticarnet.```{r}cor(datos$anticarnet, datos$edad)```Si el estudio se enfocara en la relación entre estas dos variables para analizar el efecto de sacarse el carnet antes o después sería interesante mantener las dos, pero sabiendo que en la inmensa mayoría de casos las personas se sacan el carnet alrededor de los 18 años podríamos inferir de forma bastante acertada a partir de la edad la antigüedad del carnet. Por lo tanto, eliminamos la antigüedad de carnet.```{r}datos <-subset(datos, select =-anticarnet)```Ahora podemos volver a ver el gráfico de correlaciones para comprobar que no quedan altas correlaciones:```{r}numericas <- datos[sapply(datos, is.numeric)]correlaciones <-cor(numericas, use ="pairwise.complete.obs")corrplot(correlaciones, method ="circle", type ="upper", tl.col ="black", tl.srt =45, col = colores)```Al reducir el número de variables considerando las correlaciones, obtenemos una base de datos más manejable e interpretable, lo que facilita un análisis más eficiente y una convergencia más rápida en los modelos.### 2.3 Separación de los datos en Training set y Test setPartimos al azar base en dos: 66.7% entrenamiento (667 casos) y 33.33% test (333 casos).::: {.callout-note title="Separación del dataset"}- Training set: Este conjunto se utiliza para entrenar el modelo. Contiene ejemplos históricos con características de préstamos y la información de si fueron pagados o no. El modelo analiza estos datos para identificar patrones que permitan clasificar a futuros solicitantes como buenos pagadores o posibles impagados.- Test set: Este segundo conjunto se emplea para evaluar la efectividad del modelo en datos que no ha visto antes. Sirve para comprobar si el modelo es capaz de reconocer patrones en casos nuevos y, por tanto, si sería útil para predecir con precisión en situaciones reales.:::Para separar la base de datos se establece una semilla para garantizar que los resultados sean reproducibles. Luego, los datos se ordenan pseudo-aleatoriamente, se calcula el número de muestras para el conjunto de entrenamiento y de test, y se seleccionan los datos para ambos de la database original.Después, se genera una tabla de frecuencias que muestra la distribución del número de reclamaciones en el conjunto de entrenamiento.```{r}# Permutamos los datos al azarset.seed(216514)azar <-sample(1:nrow(datos))casos <-round(nrow(datos)*0.667)train <- datos[azar[1:casos],]test <- datos[azar[(casos+1):nrow(datos)],]Tabla.Training <-table(train$numerosi)Tabla.Test <-table(test$numerosi)Tabla.Training %>%kbl(col.names =c("Número Reclamaciones", "Frecuencia"),align =rep('c', 2)) %>%kable_styling()```Realizamos la misma tabla para el test set:```{r}Tabla.Test %>%kbl(col.names =c("Número Reclamaciones", "Frecuencia"),align =rep('c', 2)) %>%kable_styling()```En ambos casos la tabla sugiere que la mayoría de las pólizas en el conjunto de datos tienen 0 o 1 reclamación. Las otras categorías tienen menos frecuencia, lo que indica que las reclamaciones mayores son poco frecuentes en esta muestra, lo cual era de esperar, e interesa mantener a la empresa.La técnica de separación usada es una de las más comunes y simples, pero puede no ser la mejor opción en todos los casos, especialmente si la distribución de las clases de los datos están desbalanceadas como acabamos de comprobar. Teniendo en cuenta esta característica de nuestra database, podemos usar el muestreo estratificado. Este método asegura que la distribución de las clases en los conjuntos de entrenamiento y test sea representativa de la distribución en el conjunto original. Es útil para evitar sesgos en la evaluación del modelo, especialmente cuando las clases menos frecuentes son importantes (que en nuestro caso lo son ya que son las que más afectan negativamente en la economía de la empresa).## 3. Análisis Exploratorio de Datos### 3.1 Variable respuestaVisualizamos la frecuencia de número de reclamaciones a lo largo del tiempo.```{r correlaciones}#| code-fold: true#| code-summary: "Mostrar código"barplot(table(train$numerosi), col = colores, border = "white", main = "Frecuencia de Número de Reclamaciones", xlab = "Número de Reclamaciones", ylab = "Frecuencia", las = 1, cex.names = 0.9, cex.axis = 0.8, cex.main = 1.2) ```Una de las propiedades más importantes de la distribución de Poisson es que su media es igual a su varianza. Ambos están representados por el parámetro λ (lambda). Esta propiedad es única y ayuda a identificar si un conjunto de datos sigue una distribución de Poisson. Esta igualdad implica que, a medida que aumenta el número esperado de sucesos, también lo hace la variabilidad en el número real de sucesos.Por lo tanto, podemos encontraar sobredispersión si la varianza de la variable de respuesta es significativamente mayor que su media, lo cual es incompatible con los supuestos del modelo Poisson como acabamos de ver.```{r}mean(train$numerosi)var(train$numerosi)```La media es 0.45, mientras que la varianza es 0.94. En un modelo Poisson ideal, la media y la varianza deberían ser iguales.La varianza es notablemente mayor que la media, lo que indica sobredispersión. Esto significa que los datos probablemente no sean bien modelados por un modelo Poisson estándar, ya que este asume que la media y la varianza son iguales. Esta igualdad implica que, a medida que aumenta el número esperado de sucesos, también lo hace la variabilidad en el número real de sucesos.El exceso de ceros ocurre cuando la proporción de ceros en la variable de respuesta es mucho mayor de lo esperado bajo una distribución Poisson.```{r}table(train$numerosi)[1] /nrow(train)```Más del 73% de las observaciones tienen un valor de cero en número de reclamaciones, lo cual esperaríamos porque indica que muchos pólizas no tienen reclamaciones, pero en términos técnicos nos puede llevar a complicaciones. Este exceso de ceros no es consistente con un modelo Poisson estándar, que generaría ceros con menor frecuencia dado el valor esperado de 0.39.Concluimos que el modelo Poisson estándar no es adecuado. Necesitarás un modelo que pueda manejar sobredispersión y exceso de ceros como por ejemplo los modelos de inflación cero. Estos utilizan una técnica estadística adaptada a los datos de recuento que tienen un exceso de resultados cero, lo cual es ideal para nuestro caso.Ahora, representamps el número de reclamaciones ajustadas según un modelo de Poisson, el número de reclamaciones que el modelo de Poisson predice que deberíamos tener para cada número de reclamaciones frente al número de reclamaciones observadas en los datos reales.```{r, warning=FALSE}#| code-fold: true#| code-summary: "Mostrar código"predict0 <- dpois(0:10, mean(train$numerosi)) * nrow(train)observado <- table(factor(train$numerosi, levels = 0:10))barplot(rbind(observado, predict0), las = 1, beside = TRUE, names.arg = c(0:9, ">=10"), col = c(colores[1], colores[5]), border = "white", main = "Distribución Poisson", xlab = "Número de Reclamaciones", ylab = "Frecuencia", las = 1, cex.names = 0.9, cex.axis = 0.8, cex.main = 1.2)```Como podemos ver, el modelo Poisson predice que deberíamos tener menos pólizas con 0 reclamaciones que las observadas, pero más de 1 reclamación, lo cual implicaría una variabilidad más constante con la media.Comparamos con una Binomial Negativa:```{r, warning=FALSE}#| code-fold: true#| code-summary: "Mostrar código"fit.bn <- goodfit(train$numerosi, type = "nbinomial", method = "MinChisq")barplot(rbind(observado, fit.bn$fitted[1:11]), las = 1, beside = TRUE, names.arg = c(0:9, ">=10"), col=c(colores[1], colores[5]), border = "white", main = "Distribución Binomial Negativa", xlab = "Número de Reclamaciones", ylab = "Frecuencia", las = 1, cex.names = 0.9, cex.axis = 0.8, cex.main = 1.2) ```Claramente la distribución Binomial Negativa se ajusta de manera más adecuada a nuestra base de datos, ya que las frecuencias ajustadas se alinean mucho mejor con las observadas, reflejando de manera más precisa la estructura de los datos.### 3.2 Relaciones variable respuesta-predictores#### 3.2.1 Tiempo de exposición al riesgoEl tiempo de exposición al riesgo se debe tratar de la siguiente manera según el contexto:- Como un offset: es el enfoque más común en modelos de conteos (Poisson, Binomial Negativa).- Como predictor: este enfoque es menos común, pero útil si creemos que el efecto del tiempo varía entre los casos.- Como peso: si queremos dar más importancia a los casos con más tiempo de exposición al riesgo.Por lo tanto, para modelos de eventos o conteos, el tratamiento más adecuado es generalmente como un offset.```{r}# Tiempo de exposición (Predictor, Peso u Offset)summary(train$duracion_poliza)```Como el tiempo de exposición afecta la tasa de siniestros, modelarlo como un offset nos permitiría ajustar por esta variable sin estimar un coeficiente independiente para ella. Al tratarlo como offset, estaríamos ajustando el modelo para que la tasa de eventos sea proporcional al tiempo de exposición, pero sin incluirla como un predictor adicional cuyo coeficiente se estimaría.Es posible que algunos contratos tengan una duración extremadamente corta, de menos de un día, como un seguro de alquiler que se contrata y se extingue el mismo día. En este caso, la "exposición al riesgo" podría ser efectivamente cero, ya que no hay un período suficiente de tiempo para que ocurra un evento asegurado (como un accidente o siniestro).```{r}sum(datos$duracion_poliza==0)sum(train$duracion_poliza==0)sum(test$duracion_poliza==0)```En nuestro caso, el conjunto de entrenamiento contiene 13 cero días al azar de los 24 registros con ceros originalmente. De la misma forma, en el conjunto de prueba hay 11 casos con duración de póliza igual a cero.En ambos conjuntos hay registros con duración de póliza de 0 días, lo cual podría ser problemático, así que los convertimos en 1, y lo repetimos para la posibilidad de que haya valores erróneos menores a 0 que se deberían a un error de imputación de los datos originales.```{r}test$duracion_poliza[test$duracion_poliza ==0] <-1test$duracion_poliza[test$duracion_poliza <0] <-1train$duracion_poliza[train$duracion_poliza ==0] <-1train$duracion_poliza[train$duracion_poliza <0] <-1```Para entender mejor la distribución de la duración de las pólizas en nuestro conjunto de datos, hemos realizado un gráfico de densidad. Este gráfico nos permite observar cómo se distribuyen los valores de duración de las pólizas, destacando posibles tendencias o patrones en la variable.```{r}#| code-fold: true#| code-summary: "Mostrar código"ggplot(train, aes(x = duracion_poliza)) +geom_density(fill = colores[1], color = colores[5], alpha =0.6) +theme_minimal() +labs(title ="Distribución de la Densidad de la Duración de la Póliza", x ="Duración de la Póliza (días)", y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) ```La distribución de la duración de las pólizas muestra un comportamiento relativamente uniforme a lo largo del eje x, con el máximo alrededor de los 175 días. Los valores mínimos se observan al inicio y hacia el final del rango, indicando una menor frecuencia de pólizas con duraciones extremadamente cortas o largas.```{r}#| code-fold: true#| code-summary: "Mostrar código"set.seed(1)uniforme <-runif(n =nrow(train), min =min(train$duracion_poliza), max =max(train$duracion_poliza))data_comparacion <-data.frame(duracion_poliza =c(train$duracion_poliza, uniforme),tipo =rep(c("Real", "Uniforme"), each =nrow(train)))ggplot(data_comparacion, aes(x = duracion_poliza, fill = tipo, color = tipo)) +geom_density(alpha =0.6) +theme_minimal() +labs(title ="Comparación de Distribución Duración Póliza y Uniforme", x ="Duración de la Póliza (días)", y ="Densidad") +scale_fill_manual(values =c(colores[1], colores2[1])) +scale_color_manual(values =c(colores[5], colores2[5])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),legend.title =element_blank(),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) +guides(fill =guide_legend(override.aes =list(alpha =0.6)))```Simulando una distribución uniforme aleatoria (con la semilla 1) podemos ver claramente que son en efecto bastante similares. Esto implica que todas las duraciones tienen una probabilidad similar de ocurrir, lo cual es consistente con una distribución uniforme, donde cada valor dentro del rango tiene la misma probabilidad.Además, puede ser interesante estudiar el límite máximo de la variable:```{r}max(train$duracion_poliza)```Que el valor máximo de la variable sea 366 días sugiere que el contrato abarca un máximo de un año calendario completo en un año bisiesto, aunque no es lo más común.```{r}cor(train$duracion_poliza, train$numerosi)```Para interpretar los valores de correlación vamos a seguir estos intervalos:::: {.callout-note title="Intervalos de interpretación del coeficiente de correlación:"}- $|r|=1$: Correlación perfecta.- $|r|\geqslant0.8$: Correlación alta.- $|r|\geqslant0.5$: Correlación media- $|r|\geqslant0.3$: Correlación débil- $|r|>0$: Correlación muy débil- $r=0$: No hay correlación.:::La correlación entre la duración de la póliza y el número de siniestros es muy débil pero positiva ($r=0.27$). Esto sugiere que, en general, a medida que aumenta la duración de la póliza, también tiende a incrementarse el número de siniestros, aunque la relación no es especialmente fuerte. Una póliza más larga brinda más tiempo durante el cual puede ocurrir un siniestro. Por tanto, el aumento en el número de siniestros es razonable desde una perspectiva temporal.De todas formas, aunque la relación existe, no es determinante. Otros factores además de la duración probablemente influyen significativamente en el número de siniestros, como el tipo de póliza, la antigüedad del vehículo o las características del asegurado. Podría ser interesante estudiar la media de siniestros para cada cantidad de días de duración.```{r}#| code-fold: true#| code-summary: "Mostrar código"df_plot <-data.frame(duracion_poliza =sort(unique(train$duracion_poliza)),promedio_siniestros =tapply(train$numerosi, train$duracion_poliza, mean))ggplot(df_plot, aes(x = duracion_poliza, y = promedio_siniestros)) +geom_line(color = colores[1], size =1.2) +geom_point(color = colores[5], size =3, shape =19) +theme_minimal() +labs(title ="Promedio de Siniestros por Duración de la Póliza",x ="Duración de la Póliza (días)", y ="Promedio de Siniestros") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"), axis.title.x =element_text(size =12),axis.title.y =element_text(size =12),axis.text =element_text(size =10), panel.grid.major =element_line(color ="gray", size =0.5), panel.grid.minor =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```Cada punto en el gráfico muestra el promedio de siniestros para una duración específica de la póliza. Parece que hay más variabilidad en el promedio de siniestros a medida que aumenta la duración de la póliza. Esto significa que, conforme pasa el tiempo, es más común que las pólizas tengan diferentes números de reclamaciones. No se observa una hay una tendencia muy clara en los puntos del gráfico, lo que sugiere que el número promedio de siniestros no tiene un patrón fácil de predecir solo con la duración de la póliza.#### 3.2.2 Potencia del vehículo```{r}summary(train$potencia)``````{r}#| code-fold: true#| code-summary: "Mostrar código"ggplot(train, aes(x = potencia)) +geom_density(fill = colores[1], color = colores[5], alpha =0.6) +theme_minimal() +labs(title ="Distribución de la Densidad de la Potencia", x ="Potencia", y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) ```Esta variable parece seguir una distribución similar a la log-normal. Lo comprobamos:```{r}#| code-fold: true#| code-summary: "Mostrar código"set.seed(1)log_normal <-rlnorm(n =nrow(train), meanlog =mean(log(train$potencia +1)), sdlog =sd(log(train$potencia +1)))data_comparacion_potencia <-data.frame(potencia =c(train$potencia, log_normal),tipo =rep(c("Real", "Log-Normal"), each =nrow(train)))ggplot() +geom_density(data =subset(data_comparacion_potencia, tipo =="Real"), aes(x = potencia, fill ="Real"), alpha =0.6, color = colores[5]) +geom_density(data =subset(data_comparacion_potencia, tipo =="Log-Normal"), aes(x = potencia, fill ="Log-Normal"), alpha =0.4, color = colores2[5]) +theme_minimal() +labs(title ="Comparación de Distribución Potencia y Log-Normal", x ="Potencia", y ="Densidad") +scale_fill_manual(values =c("Real"= colores[1], "Log-Normal"= colores2[1])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),legend.title =element_blank(),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) +guides(fill =guide_legend(override.aes =list(alpha =0.6)))```La distribución real de la potencia y la simulada con una log-normal son muy similares, lo que sugiere que la potencia de los vehículos en el conjunto de datos podría aproximarse razonablemente a este tipo de distribución. Esto puede darse porque en muchas situaciones reales, las variables relacionadas con características físicas o de rendimiento (como la potencia del motor) tienden a distribuirse de forma asimétrica, con una mayoría de valores en un rango medio y una cola hacia valores más altos que en este caso serían los coches de carreras ode lujo. Este comportamiento es característico de una distribución log-normal.Ahora estudiamoscómo varía el número promedio de siniestros (variable respuesta) en función de la potencia del vehículo:```{r}medias <-tapply(train$numerosi, train$potencia, mean)``````{r}#| code-fold: true#| code-summary: "Mostrar código"plot(sort(unique(train$potencia)), medias, pch =20, col = colores[3], xlab ="Potencia del vehículo", ylab ="Número promedio de siniestros", main ="Relación entre Potencia y Número Promedio de Siniestros", cex =1.2, # Tamaño de los puntoscex.axis =1, # Tamaño de los ejescex.lab =1.1, # Tamaño de las etiquetascex.main =1.3, # Tamaño del títulolwd =2, # Grosor de la línea de los puntoslty =1, # Tipo de línea (sin línea, solo puntos)xlim =c(min(sort(unique(train$potencia))) -10, max(sort(unique(train$potencia))) +10))```Aparentemente, la mayoría de los puntos están agrupados en los valores más bajos de potencia y número promedio de siniestros, lo que sugiere que los vehículos con menor potencia tienden a tener menos siniestros. A parte, hay algunos puntos fuera de la masa principal de puntos, con valores más altos tanto en potencia como en el número promedio de siniestros. Esto indica que, aunque no es la norma, algunos vehículos con alta potencia tienen un número significativo de siniestros.En nuestra base de datos no parece haber una tendencia específica, ya que los puntos se distribuyen al rededor de 0. Esto podría sugerir que los vehículos de menor potencia son más seguros o que son conducidos de manera más conservadora.```{r}cor(train$potencia, train$numerosi)```Volviendo a la interpretación de cada rango de correlación, un valor de $r=0.08$ indica una relación lineal directa muy débil, lo cual podíamos esperar del gráfico anterior ya que no se observa ninguna relación clara.La discretización de una variable continua, como en este caso la potencia de un vehículo, implica dividirla en grupos o categorías. Este proceso puede ser útil para facilitar la interpretación de la relación entre esa variable y otras, como el número de siniestros. Se pueden usar los cuantiles de la variable para dividir los datos en grupos más equitativos en términos de cantidad de observaciones.```{r}# Discretizar la potencia usando cuantilesqnt <-quantile(train$potencia, seq(0, 1, .2))qnt[1] <- qnt[1] -0.01# Ajustar el primer cuantílqnt[5] <- qnt[5] +0.01# Ajustar el quinto cuantíltrain$g.potencia <-cut(train$potencia, breaks = qnt)train$g.potencia <-factor(train$g.potencia, levels =levels(train$g.potencia))medias <-tapply(train$numerosi, train$g.potencia, mean)``````{r}#| code-fold: true#| code-summary: "Mostrar código"ggplot(data =data.frame(x =1:length(medias), y = medias), aes(x = x, y = y)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5) +scale_x_continuous(breaks =1:length(medias), labels =levels(train$g.potencia)) +labs(title ="Promedio de Siniestros por Grupo de Potencia", x ="Rango de Potencia", y ="Promedio de Siniestros") +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```Como podemos ver, los vehículos en el rango de potencia (90-112) parecen estar asociados con un mayor número de siniestros. Esto podría deberse a que estos vehículos tienen suficiente potencia para ser utilizados de manera más agresiva o en condiciones que aumentan el riesgo de siniestros.Los vehículos en el rango de potencia (35-90) y (112-272) tienen menos siniestros en promedio, lo cual podría sugerir que estos vehículos son más seguros o que son conducidos de manera más conservadora, de hecho probablemente estos vehículos son tractores, patinetes y otros similares.La variabilidad en otros rangos muestra que hay diferencias aparentemente significativas en el número promedio de siniestros según la potencia del vehículo.La discretización de la variable ha sido muy interesante, podemos estar interesados en hacerlo de nuevo, por eso hemos creado una función para hacerlo de forma más automática:```{r}#| code-fold: true#| code-summary: "Mostrar código"plot.dist <-function(predictor, cortes=8, variable=F){ qnt <-quantile(predictor, seq(0,1,1/cortes)) qnt[1] <- qnt[1]-0.01; qnt[cortes] <- qnt[cortes]+0.01 V.discreta <-cut (predictor, breaks= qnt) medias <-tapply(train$numerosi, V.discreta, mean)if (variable == T) return(V.discreta)} # Fin función```Esta función nos ayudará a realizar los cálculos de discretización de forma mucho más rápida.#### 3.2.3 Antigüedad del vehículoEstudiamos ahora la variable de antigüedad del vehículo.```{r}#| code-fold: true#| code-summary: "Mostrar código"ggplot(train, aes(x = antivehi)) +geom_density(bw =0.5, fill = colores[3], color = colores[5], alpha =0.6) +labs(title ="Distribución de la Antigüedad del Vehículo",x ="Antigüedad del Vehículo (años)",y ="Densidad" ) +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12),axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```El gráfico no parece ajustarse a una distribución específica, pero sí muestra una clara tendencia decreciente. Observamos que hay muchos vehículos relativamente nuevos, con menos de 5 años de antigüedad, seguido de una disminución significativa a partir de entre los 11 y 12 años.Para realizar la comparación, identificamos que la distribución más similar a nuestros datos es la chi-cuadrado, aunque no parece tan similar como en los casos anteriores. Esta distribución se caracteriza por su asimetría, comenzando con valores altos que disminuyen progresivamente. Es adecuada en este caso porque refleja una tendencia inicial elevada que va decreciendo, algo similar a lo observado en nuestros datos de antigüedad.```{r}#| code-fold: true#| code-summary: "Mostrar código"scale_chi <-mean(train$antivehi) /2set.seed(1)chi_data <-rchisq(n =length(train$antivehi), df =2) * scale_chidata_comparacion_antivehi <-data.frame(antiguedad =c(train$antivehi, chi_data),tipo =c(rep("Real", length(train$antivehi)), rep("Chi-Cuadrado", length(chi_data))))ggplot() +geom_density(data =subset(data_comparacion_antivehi, tipo =="Real"), aes(x = antiguedad, fill ="Real"), alpha =0.6, color = colores[5]) +geom_density(data =subset(data_comparacion_antivehi, tipo =="Chi-Cuadrado"), aes(x = antiguedad, fill ="Chi-Cuadrado"), alpha =0.4, color = colores2[5]) +theme_minimal() +labs(title ="Comparación de Distribución Antigüedad del Vehículo y Chi-Cuadrado", x ="Antigüedad del Vehículo (años)", y ="Densidad") +scale_fill_manual(values =c("Real"= colores[1], "Chi-Cuadrado"= colores2[1])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),legend.title =element_blank(),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm")) +guides(fill =guide_legend(override.aes =list(alpha =0.6)))```Como podemos observar, ambas distribuciones son parecidas, salvo por un repunte en los datos reales alrededor de los 10 años de antigüedad. Este comportamiento podría estar relacionado con un mayor número de vehículos que alcanzan esa antigüedad debido a características del mercado o patrones de uso. A pesar de esta diferencia, la similitud se recupera al volver a descender posteriormente.```{r}#| code-fold: true#| code-summary: "Mostrar código"antivehi_data <-data.frame(antivehi =sort(unique(train$antivehi)),mean_numerosi =tapply(train$numerosi, train$antivehi, mean))ggplot(antivehi_data, aes(x = antivehi, y = mean_numerosi)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5) +theme_minimal() +labs(title ="Relación entre Antigüedad del Vehículo y Número de Siniestros",x ="Antigüedad del Vehículo (años)",y ="Promedio de Número de Siniestros" ) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```Los vehículos más nuevos (0 años) tienen el promedio de siniestros más alto, más de 0.87. Esto podría deberse a la falta de experiencia del conductor con el vehículo nuevo. A medida que los vehículos son ligeramente más antiguos (de 1 a 5 años), el promedio de siniestros disminuye, alcanzando un mínimo local a los 7 años. Después, el promedio de siniestros muestra algunas fluctuaciones, pero en general vemos que decrece con algunos incrementos ocasionales.Tras ver el gráfico, puede ser interesante ver la correlación entre la variable de respuesta, número de siniestros, y la antigüedad del vehículo.```{r}cor(train$antivehi, train$numerosi) ```Hay una muy débil correlación directa ($r=-0.18$), es decir, existe una pequeña tendencia a que, a medida que aumenta la antigüedad del vehículo, disminuye el número de siniestros en promedio. El hecho de que sea tan débil implica que la antigüedad del vehículo no es un predictor particularmente fuerte del número de siniestros.#### 3.2.4 Tipo de vehículoLa variable tipo de vehículo es de tipo factor, por lo tanto es preciso otro tipo de análisis.```{r}summary(train$tipovehi)```Estudiamos su estructura con un gráfico:```{r}#| code-fold: true#| code-summary: "Mostrar código"ggplot(train, aes(x = tipovehi, y = numerosi)) +stat_summary(fun = mean, geom ="bar", fill = colores[1], color = colores[5], width =0.6) +theme_minimal() +labs(title ="Promedio de Siniestros por Tipo de Vehículo", x ="Tipo de Vehículo", y ="Promedio de Siniestros") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),axis.text.x =element_text(angle =45, hjust =1), panel.border =element_blank(),plot.margin =margin(2, 1, 1, 1, "cm"))```El gráfico muestra el promedio de siniestros (reclamaciones) para diferentes tipos de vehículos. Los vehículos del tipo 120 y 310 tienen los promedios de siniestros más altos (aproximadamente 1 y 0.65 respectivamente). Esto sugiere que estos tipos de vehículos tienen más probabilidades de estar involucrados en siniestros. Los vehículos del tipo 250 y 300 no reportan a penas siniestros y los del resto de tipos tienen un promedio bajo (menos de 0.5). Esto indica que estos vehículos son menos propensos a involucrarse en siniestros.Hay diferencias significativas en el promedio de siniestros entre los diferentes tipos de vehículos, lo cual es importante para los estudios de evaluación de riesgos en seguros. Estas diferencias hacen que la distribución no se asemeje a ninguna distribución común.Ahora, al ser la variable de tipo factor, comprobamos si hay diferencias significativas en el número de siniestros entre los distintos tipos de vehículos. Usando ANOVA sabemos si tipo de vehículo tiene un impacto importante en el número de siniestros o si es más probable que las diferencias que vemos sean solo por casualidad.```{r}mod1 <-lm(train$numerosi ~ train$tipovehi)anova(mod1)```::: {.callout-note title="Umbrales comunes de significancia"}El coeficiente de correlación mide la intensidad y dirección de la relación lineal entre dos variables:- $p-value=0.05$: Este es el umbral más comúnmente utilizado.- $p-value=0.01$: Un umbral más estricto.- $p-value=0.1$: A veces se usa en investigaciones exploratorias, pero no se considera concluyente.:::El análisis de varianza muestra que el tipo de vehículo no tiene un impacto significativo en el número de siniestros. El valor F es bajo y el valor p es alto (0.5675), no entra en ninguno de los casos mencionados, lo que sugiere que las diferencias en el número de siniestros entre los diferentes tipos de vehículos podrían ser atribuibles al azar y no son estadísticamente significativas. En otras palabras, el tipo de vehículo no parece ser un factor importante para predecir el número de siniestros en este caso.```{r}confint(mod1)```Aquí vemos los intervalos de confianza al 95% para los coeficientes estimados del modelo lineal que hemos creado. Específicamente, estos números nos dicen el rango dentro del cual creemos que se encuentra el valor real del coeficiente, con un 95% de certeza.Para todos los tipos de vehículo mencionados, los intervalos de confianza incluyen el cero en su rango, lo que indica que, en efecto, no podemos afirmar que exista una relación significativa entre estos tipos de vehículos y el número de siniestros. En otras palabras, no podemos concluir que el tipo de vehículo tenga un impacto claro sobre el número de siniestros en el modelo.#### 3.2.5 Plazas de vehículoEstudiamos ahora la variable plazas de vehículo.```{r}table(train$plazasvehi)```Como vemos, es una variable muy desbalanceada, aunque en términos generales, plazasvehi tiene dos categorías (5 y 6 plazas), la gran mayoría de los vehículos tiene 5 plazas, mientras que solo una pequeña proporción tiene 6 plazas.```{r}tapply(train$numerosi, train$plazasvehi, mean)```Esta funciónnos permite ver cómo varía el promedio de siniestros (numerosi) según el número de plazas del vehículo (plazasvehi), lo cual elimina el efecto del desbalance. El resultado nos indica que para los vehículos con 5 plazas, el promedio de los siniestros (train\$numerosi) es 0.4471, y para los vehículos con 6 plazas, el promedio de los siniestros es 0.6.Los vehículos con 6 plazas suelen ser más grandes, lo que podría implicar que son utilizados de manera diferente o para viajes más largos, lo que podría aumentar las probabilidades de que ocurran siniestros. En cambio, los vehículos más pequeños de 5 plazas a menudo son conducidos de manera diferente o menos veces, lo que podría reducir la exposición al riesgo de siniestros.#### 3.2.6 EdadAhora estudiamos la variable de edad. Realizamos de nuevo el gráfico de la distribución de la variable:```{r}ggplot(train, aes(x = edad)) +geom_density(fill = colores[3], color = colores[5], alpha =0.6, bw =0.5) +theme_minimal() +labs(title ="Distribución de la Edad", x ="Edad", y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```La distribución de la edad sigue una forma que se asemeja bastante a una distribución normal, lo cual es común en los conjuntos de datos con esta variable. Las vemos constrastadas en el siguiente gráfico:```{r}#| code-fold: true#| code-summary: "Mostrar código"media_edad <-mean(train$edad, na.rm =TRUE)sd_edad <-sd(train$edad, na.rm =TRUE)ggplot(train, aes(x = edad)) +geom_density(aes(fill ="Datos Reales"), color = colores[5], alpha =0.6, bw =0.5) +stat_function(fun = dnorm, args =list(mean = media_edad, sd = sd_edad), color = colores2[5], size =0.6,geom ="area", aes(fill ="Normal Ajustada"), alpha =0.3) +theme_minimal() +labs(title ="Comparación de Distribución Normal y Edad", x ="Edad", y ="Densidad") +scale_fill_manual(name ="Distribuciones", values =c("Datos Reales"= colores[3], "Normal Ajustada"= colores2[3])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```La distribución observada en los datos reales sigue una forma similar a una distribución normal, pero presenta algunas irregularidades. El gráfico muestra un pico notable cerca de los 50 años, lo que indica que la mayoría de los casos se concentran en este rango de edad. Además, se observa un descenso en la densidad justo después de los 50 años, lo que sugiere una disminución en la frecuencia de los datos en ese grupo, para el cual no queda muy clara la interpretación.Estudiamos ahora su relación con la variable respuesta:```{r}cor(train$edad, train$numerosi) ```Para empezar, el coeficiente de correlación (\$r=-0.008) nos indica una correlación muy débil negativa. A medida que la edad del conductor aumenta, el número de siniestros disminuye en sentido general, pero es tan débil la relación que no nos da ninguna información relevante.Ahora vemos esta relación graficada:```{r}#| code-fold: true#| code-summary: "Mostrar código"edad_data <-data.frame(edad =sort(unique(train$edad)),mean_numerosi =tapply(train$numerosi, train$edad, mean))ggplot(edad_data, aes(x = edad, y = mean_numerosi)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5) +theme_minimal() +labs(title ="Relación entre Edad y Número de Siniestros",x ="Edad",y ="Promedio de Número de Siniestros" ) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```En efecto, no queda muy clara la relación entre las variables. Por lo tanto, tal como hemos aprendido, discretizar la variable puede ayudar a su interpretación. La discretizamos usando la función que habíamos creado previamente.```{r}dist.test <-function(predictor.train, predictor.test, cortes=10) { qnt <-quantile(predictor.train, seq(0, 1, 1/cortes)) qnt[1] <- qnt[1] -0.01 qnt[cortes] <- qnt[cortes] +0.01 V.discreta <-cut(predictor.test, breaks = qnt)return(V.discreta)}train$edad_discretizada <-dist.test(train$edad, train$edad, cortes =10)test$edad_discretizada <-dist.test(train$edad, test$edad, cortes =10)levels(test$edad_discretizada) <-levels(train$edad_discretizada)``````{r}#| code-fold: true#| code-summary: "Mostrar código"edad_data_discretizada <-data.frame(edad_discretizada =levels(train$edad_discretizada),mean_numerosi =tapply(train$numerosi, train$edad_discretizada, mean))edad_data_discretizada$edad_discretizada <-factor(edad_data_discretizada$edad_discretizada, levels = edad_data_discretizada$edad_discretizada)ggplot(edad_data_discretizada, aes(x =as.numeric(edad_discretizada), y = mean_numerosi)) +geom_point(color = colores[5], size =2) +geom_line(color = colores[3], size =0.5, aes(group =1)) +scale_x_continuous(breaks =1:length(edad_data_discretizada$edad_discretizada), labels = edad_data_discretizada$edad_discretizada) +theme_minimal() +labs(title ="Relación entre Edad Discretizada y Promedio de Número de Siniestros",x ="Grupo de Edad Discretizado",y ="Promedio de Número de Siniestros" ) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```Ahora la relación queda mucho más clara. Los conductores más jóvenes (36-40 años) tienen un promedio de siniestros bastante alto, probablemente debido a la falta de experiencia al volante y comportamientos de conducción más arriesgados. Los conductores de edad joven-media (20-36) muestra una disminución en el promedio de número de siniestros, lo cual indica un mejor manejo al volante o hábitos más seguros al volante con los años. En los conductores de edad media (40-53) disminuye el promedio de siniestros, lo cual puede deberse a que los conductores en este grupo de edad, además de tener más experiencia al volante, a menudo tienen una vida profesional y familiar activa, por lo que podría significar que conducen con más cuidado. Los conductores más mayores (53-84) tienen un promedio de siniestros más bajo, posiblemente porque conducen con mayor precaución o debido a que no suelen conducir largos y complicados viajes sino cortos viajes del día a día para lo esencial.#### 3.2.7 Sexo conductorEn este caso analizamos la variable del sexo del conductor. Esta variable vuelve a ser de tipo factor así que la vamos a analizar de forma parecida a la variable de tipo ed vehículo.```{r}table(train$sexocondu)```De acuerdo con los datos, la mayoría de las conductoras son femeninas (dando por hecho que M es masculino y V femenino), ya que hay 536 observaciones de conductores femeninos en comparación con 131 de conductores masculinos. Esto indica que el grupo de conductores en este conjunto de datos está desbalanceado en favor de las mujeres.```{r}tapply(train$numerosi, train$sexocondu, mean)```Como vemos, los hombres en esta base de datos tienden a tener más siniestros en comparación con las mujeres, pero debemos comprobar si hay diferencias significativas en el número de siniestros entre estos. Usando ANOVA sabemos si el sexo tiene un impacto importante en el número de siniestros o si es más probable que las diferencias que vemos sean solo por casualidad.```{r}mod2 <-lm(train$numerosi ~ train$sexocondu)anova(mod2)```Recordando los umbrales comunes de significancia, el análisis de varianza muestra que el sexo tiene un impacto significativo al nivel de significatividad 0.01 en el número de siniestros. El valor p es bajo (0.1244), lo que sugiere que las diferencias en el número de siniestros entre los dos sexos podrían deberse. La variable sexo del conductor parece ser un factor importante para predecir el número de siniestros.```{r}confint(mod2)```Este intervalo de confianza sugiere que con un 95% de confianza, el efecto de ser mujer en la variable objetivo es negativo. Es decir, el ser mujer reduce el número de siniestros en comparación con los hombres tal como esperábamos con los estudios anteriores.#### 3.2.8 Antigüedad en la compañíaAhora estudiamos la variable antigüedad en la compañía:```{r}#| code-fold: true#| code-summary: "Mostrar código"ggplot(train, aes(x = anticia)) +geom_density(fill = colores[3], color = colores[5], alpha =0.6, bw =0.5) +theme_minimal() +labs(title ="Distribución de la Antigüedad en la Compañía",x ="Antigüedad del Vehículo (años)",y ="Densidad") +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```La mayoría de los clientes en la muestra llevan entre 0 y 5 años en la compañía Este rango tiene el pico de densidad más alto, cerca de 0.20, lo que sugiere que los clientes nuevos son los más comunes. A medida que la antigüedad aumenta, la densidad disminuye gradualmente. Esta distribución nos recuerda a una de tipo exponencial. Las vemos comparadas a continuación:```{r}#| code-fold: true#| code-summary: "Mostrar código"lambda_anticia <-1/mean(train$anticia, na.rm =TRUE)ggplot(train, aes(x = anticia)) +geom_density(aes(fill ="Datos Reales"), color = colores[5], alpha =0.6, bw =0.5) +# Densidad observadastat_function(fun = dexp, args =list(rate = lambda_anticia), color = colores2[5], size =0.6,geom ="area", aes(fill ="Exponencial Ajustada"), alpha =0.3) +# Curva exponencial ajustadatheme_minimal() +labs(title ="Comparación Distribución Exponencial y Antigüedad del Vehículo", x ="Antigüedad del Vehículo (años)", y ="Densidad") +scale_fill_manual(name ="Distribuciones", values =c("Datos Reales"= colores[3], "Exponencial Ajustada"= colores2[3])) +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```Como podemos observar, las distribuciones son bastante similares. Sin embargo, la distribución real no comienza tan alta como la exponencial y presenta un descenso más pronunciado entre aproximadamente los 3 y los 7 años. Además, hacia el final, la distribución real no muestra una caída tan pronunciada, sino más mantenida. Aunque las dos distribuciones son similares, las diferencias en el inicio, el comportamiento intermedio y el final indican que la distribución real tiene características particulares que no se ajustan perfectamente al modelo exponencial, sugiriendo que otros factores están influyendo en la antigüedad de los clientes en la compañía.Ahora estudiamos su relación con la variable respuesta:```{r}cor(train$anticia, train$numerosi)```Observamos una correlación muy débil ($r=-0.01$), lo que sugiere que, en general, a medida que un cliente lleva más años con la empresa, tiende a tener menos accidentes con el vehículo. Sin embargo, esta relación es tan débil que no es lo suficientemente "significativa" como para considerarse relevante.```{r}#| code-fold: true#| code-summary: "Mostrar código"anticia_data <-data.frame(anticia =sort(unique(train$anticia)),mean_numerosi =tapply(train$numerosi, train$anticia, mean))ggplot(anticia_data, aes(x = anticia, y = mean_numerosi)) +geom_point(color = colores[5], size =2) +# Puntos en el gráficogeom_line(color = colores[3], size =0.7, aes(group =1)) +# Línea de conexión entre puntoslabs(title ="Relación entre Antigüedad del Vehículo y Promedio de Siniestros",x ="Antigüedad del Vehículo (años)",y ="Promedio de Número de Siniestros" ) +theme_minimal() +# Estilo minimalistatheme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),panel.border =element_blank(),plot.margin =margin(1, 1, 1, 1, "cm") )```De nuevo, el gráfico no deja clara ninguna relación entre las variables, así que discretizamos la variable agrupándola en si lleva menos de 3 años en la empresa (nuevo o new) o más (veterano o vet).```{r}train$g.anticia <-ifelse(train$anticia <3, "New", "Vet") train$g.anticia <-as.factor(train$g.anticia)tapply(train$numerosi, train$g.anticia, mean)```A pesar de que los vehículos nuevos tienen un promedio ligeramente mayor de siniestros que los vehículos veteranos, la diferencia es bastante pequeña. Esto sugiere que la antigüedad del vehículo, según esta discretización en "New" y "Vet", no parece tener un impacto significativo en el número promedio de siniestros.#### 3.2.9 Usos del VehículoPor último, analizamos la variable uso del vehículo que ya habíamos discretizado en dos grupos al inicio de la práctica.```{r}table(train$usovehi2)tapply(train$numerosi, train$usovehi2, mean)```Los vehículos con un uso más alto (grupo G1 con uso mayor a 500) tienen un promedio de siniestros más alto en comparación con los vehículos con un uso menor (grupo G2 con uso menor a 100). Esto podría sugerir que los vehículos con más uso (G1) están involucrados en más accidentes o siniestros en promedio.La convertimos a tipo factor para hacer un análisis más profundo de estas diferencias:```{r}train$usovehi2 <-as.factor(train$usovehi2)mod3 <-lm(train$numerosi ~ train$usovehi2)anova(mod3)```Dado que el valor p es alto, 0.6221, no hay evidencia estadística suficiente para concluir que el uso del vehículo tiene un efecto significativo sobre el número de siniestros. En otras palabras, las diferencias entre los grupos G1 y G2 en términos de su promedio de siniestros podrían darse simplemente por el azar.```{r}confint(mod3)```El intervalo de confianza que incluye el valor 0 indica que la diferencia en el efecto entre el grupo G1 y el grupo G2 no es estadísticamente significativa. Por lo tanto, como ya hemos visto, no podemos afirmar que ser parte del grupo G2 tenga un impacto significativo en el número de siniestros.### 3.3 Algunas conclusiones del EDAEn el análisis exploratorio de datos (EDA), la mayoría de las variables estudiadas y sus respectivas transformaciones parecen tener una capacidad predictiva notable por sí solas, con algunas excepciones como las variables plazas del vehículo y tipo de vehículo, que no muestran un impacto significativo. Por tanto, se optará por incluir todas estas variables en el modelo inicial.En el caso de edad o antigüedad del carnet, hemos observado una relación no lineal con la variable dependiente. Este comportamiento sugiere que la edad podría no seguir una tendencia lineal simple, y es por eso que hemos decidido incluir transformaciones de la edad, como edad al cuadrado y edad al cubo, que pueden captar esta no linealidad de manera más efectiva, que funciona como alternativa al uso de los modelos GAM.```{r}# Edad cuadrado y cubotrain$edad2 <- train$edad^2train$edad3 <- train$edad^3```Hay muchas oportunidades para construir más predictores que podrían mejorar el modelo. Algunas posibles ideas incluyen:- Discretizar variables adicionales o hacerlo de manera diferente, por ejemplo, agrupando rangos de edad, antigüedad del vehículo, etc.- Incorporar efectos de interacción entre variables, ya que algunas combinaciones de variables podrían tener un impacto mayor al combinarse.- Crear nuevas variables derivadas a partir de las que ya tenemos disponibles. Un ejemplo sería utilizar el código postal para generar una nueva variable que capture la localización geográfica y su influencia potencial en el comportamiento de los siniestros o cualquier otra variable de interés.Esto puede enriquecer el modelo y mejorar su capacidad predictiva, por lo que seguir explorando nuevas combinaciones y transformaciones de variables sería una parte clave del proceso de modelado en otro caso.## 4 Evaluando modelos**Propuestas de modelización**Especificación:- Variables que van a entrar en el predictor lineal- Distribución para la variable respuesta- Función link- Y con enfoque ML, tipo de regularización/penalización**Distribuciones para la variable respuesta**Como modelos para la variable respuesta podemos testar:- Poisson- Binomial Negativa- ZI (inflada en cero) + Poisson- ZI (inflada en cero) + BN- Normal **Funciones link**Como funciones link:- Logaritmo<br>- Identidad **Regularización** Penalización: \* Net Elastic (Lasso, Ridge)::: {.callout-note title="Modelos sin y con regularización"}- Modelos sin regularización: Estos modelos ajustan los datos de forma directa, maximizando el ajuste sin restricciones. Aunque pueden ofrecer un buen desempeño inicial, son más propensos al sobreajuste, capturando ruido en lugar de patrones generales.- Modelos con regularización: Incorporan penalizaciones que limitan la complejidad del modelo, favoreciendo soluciones más simples y generalizables. Ayudan a evitar el sobreajuste y suelen ser más robustos en escenarios reales..:::### 4.1 Modelos sin regularización Cuando utilizamos modelos sin regularización hemos de estar atentos a los problemas de multicolinealidad.Antes de comenzar a modelizar hemos de eliminar algunas variables por que es lo razonable, lógico, desde el punto de vista de modelización y de computación.Eliminamos codigopostal y ccaa ya que no aportan información relevante al modelo y no necesariamente tienen una relación significativa o directa con el número de siniestros, ya que ambas variable son identificadores geográficos y podrían afectar a la complejidad del modelo en lugar de mejorarlo. También eliminamos las variables iniciop y finalp, debido a que ya hemos calculado con anterioridad la diferencia entre estas para determinar la duración de la póliza, por lo que mantener estas variables sería redundante y podría dificultar la interpretación del modelo.```{r}train <-subset(train, select =-c(codigopostal, ccaa, iniciop, finalp))```**Ajuste modelos (clásicos) de Poisson**```{r}Po.log <-glm(numerosi ~ . , data = train, family =poisson(link ="log"))# Obtenemos los valores iniciales y nos aseguramos que no haya coeficientes negativos o NAstarting.values <-coef(Po.log)starting.values[starting.values <0] <-0.0001starting.values[is.na(starting.values)] <-0.0001# Ajuste modelo Poisson con enlace identidad usando los valores inicialesPo.id <-glm(numerosi ~ . , data = train, start = starting.values,family =poisson(link ="identity"))```**Ajuste modelo (clásico) Binomial Negativo**```{r}BN.log <-glm.nb(numerosi ~ . , data = train, link = log)# starting.values <- coef(BN.log)# starting.values[starting.values < 0] <- 0.0001# BN.id <- glm.nb(numerosi ~ . , data = train,# start = starting.values, link = identity)```**Ajustes de modelos (clásicos) inflados de ceros**```{r}# Modelo inflado de ceros con PoissonZI.Po <-zeroinfl(numerosi ~ ., data = train, dist ="poisson")# Modelo inflado de ceros con Binomial NegativaZI.BN <-zeroinfl(numerosi ~ ., data = train, dist ="negbin")```**Normal**La variable de respuesta, numerosi, representa el número de siniestros para cada observación, por lo que se trata de una variable discreta que toma valores enteros no negativos entre los cuales es común observar una gran cantidad de ceros.Debido a esto, realizar un ajuste con modelo Normal no tiene sentido ya que la Normal es una distribución que admite valores negativos y podría generar errores significativos en la estimación de los parámetros y predicciones poco precisas. Por lo tanto, descartamos este modelo.#### 4.1.1 El rol del tiempo de exposición: Modelos con Offsets (y pesos)A continuación, realizamos diversos modelos para analizar el número de siniestros, incluyendo el manejo del tiempo de exposición utilizando offsets y pesos.```{r}#| code-fold: true#| code-summary: "Mostrar código"# Poisson con offsetPo.log.off <-glm(numerosi ~ . ,data =subset(train, select =-c(duracion_poliza)),offset =log(train$duracion_poliza /365.25), # Convertimos la duración a años family =poisson(link ="log"))# Poisson con pesos (corrigiendo `weight` a `weights`)Po.log.wgt <-glm(numerosi ~ . ,data =subset(train, select =-c(duracion_poliza)),weights = train$duracion_poliza /365.25, family =poisson(link ="log"))```El término offset nos permite incorporar la duración de la póliza como una variable que ajusta la tasa de siniestros por unidad de tiempo, garantizando que la duración de la póliza se considere de manera proporcional al análisis. En este caso, la duración se convierte a años, y su efecto se incluye de manera multiplicativa en el modelo de Poisson, ajustando las tasas de siniestros en función del tiempo de exposición al riesgo. Por otro lado, el modelo con pesos utiliza la duración de la póliza como un peso, asignando mayor relevancia a las observaciones con una mayor duración, y corrigiendo así el modelo para que tenga en cuenta de manera adecuada la variable de tiempo sin necesidad de usar un término offset explícitamente.```{r}# Comparación de AIC'sAIC(Po.log.off)AIC(Po.log.wgt) *nrow(train) /sum(train$duracion_poliza /365.25)```Observamos que el modelo con offset presenta un AIC más bajo que el modelo con pesos. Esto sugiere que el primer modelo se ajusta mejor a los datos, además de ser más interpretable en términos de tasas de siniestros ajustadas por la duración de la póliza.```{r}#| code-fold: true#| code-summary: "Mostrar código"# Poisson con offset (especificación alternativa)of_var <-log(train$duracion_poliza /365.25)formula <- numerosi ~ potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + duracion_poliza + usovehi2 + g.potencia +offset(of_var) + g.anticia + edad2 + edad3 + edad_discretizadaPo.log.off2 <-glm(formula, data = train, family = poisson)# Binomial Negativa con offsetBN.log.off <-glm.nb(formula, data = train, link = log)# Modelos inflados de ceros, con offsetformula2 <-formula(numerosi ~ potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia +offset(log(duracion_poliza/365.25)) + g.anticia + edad2 + edad3| potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia + g.anticia + edad2 + edad3 + edad_discretizada)ZI.Po.off <-zeroinfl(formula2, data = train, dist ="poisson")ZI.BN.off <-zeroinfl(formula2, data = train, dist ="negbin")```Realizamos otros modelos utilizando una fórmula más detallada, especificando los predictores. La variable respuesta aparece como un offset explícito, lo cual aclara su contribución al modelo.La Binomial Negativa permite manejar la sobredispersión que puede llegar a ocurrir con Poisson, y los modelos inflados de cero se utilizan cuando hay una gran cantidad de ceros en la variable dependiente, en este caso el número de siniestros, condición que no puede ser capturada adecuadamente por Poisson o Binomial Negativa.#### 4.1.2 Evaluación clásica de los modelosProcedemos a comparar los diferentes modelos predictivos utilizando el criterio AIC.```{r}AICs <-c(AIC(Po.log), AIC(Po.id), AIC(BN.log), AIC(ZI.Po), AIC(ZI.BN),AIC(Po.log.off), # ‘Ajuste’ por número de observacionesAIC(Po.log.wgt)*nrow(train)/sum(train$duracion_poliza/365.25),AIC(BN.log.off), AIC(ZI.Po.off), AIC(ZI.BN.off))``````{r}#| code-fold: true#| code-summary: "Mostrar código"errores_clasicos <-data.frame(modelo =c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", "Po.log.off", "Po.log.wgt", "BN.log.off", "ZI.Po.off", "ZI.BN.off"), AICs)errores_clasicos %>%kbl(align =rep('c', 3)) %>%kable_styling()``````{r}#| code-fold: true#| code-summary: "Mostrar código"colores1 <-carto_pal(11, "Teal")barplot( errores_clasicos$AICs, names.arg = errores_clasicos$modelo, col = colores1[rank(errores_clasicos$AICs)], border ="white", main ="Valores AIC por Modelo", xlab ="Modelos", ylab ="AIC", las =1, cex.names =0.6,cex.axis =0.8, cex.main =1.2)```Vemos que el modelo Poisson inflado de ceros con offset presenta el AIC más bajo (AIC = 1098.081), lo cual nos indica que es el modelo que mejor se ajusta a los datos. Elegiríamos este modelo, ya que un AIC bajo sugiere que este tiene una mejor relación entre la capacidad de ajuste y el número de parámetros, evitando así el sobreajuste.#### 4.1.3 Evaluación enfoque MLDefinición de nuevos predictores definidos en el conjunto de training en test. ```{r}#| code-fold: true#| code-summary: "Mostrar código"qnt <-quantile(train$potencia,seq(0,1,.2))qnt[1] <- qnt[1]-0.01qnt[5] <- qnt[5]+0.01test$g.potencia <-cut(test$potencia, breaks= qnt)test$g.anticia <-ifelse(test$anticia <3, "New", "Vet")test$edad2 <- test$edad^2test$edad3 <- test$edad^3test <-subset(test, select =-c(codigopostal, ccaa, iniciop, finalp))test <- test[test$duracion_poliza !=0,]```Discretizamos la potencia usando cuantiles en el test, como ya hicimos anteriormente en el training set. También clasificamos la antigüedad de la compañía en "Nueva" o "Veterana", creamos tranformaciones para capturar relaciones no lineales entre edad y otras variables en modelos predictivos, eliminamos columnas innecesarias y por último, filtramos los datos para eliminar observaciones donde la duración de la póliza sea 0 para evitar errores posteriores. ::: {.callout-warning title="Eliminación de observaciones"}Tras haber intentado generar las predicciones, obtuvimos un error sobre los niveles 5 y 19 de la variable prov. En esta tabla observamos que dichos niveles no están presentes en el conjunto de entrenamiento mientras que sí lo están en el conjunto de prueba, generando un error ya que los modelos no sabrán tratar estos niveles. Para evitar errores durante la generación de predicciones y demás, eliminamos las observaciones de dichos niveles en el test set.:::```{r}#| code-fold: true#| code-summary: "Mostrar código"resumen <-data.frame(Prov =names(table(train$prov)),Train =as.integer(table(train$prov)),Test =as.integer(table(test$prov)[names(table(train$prov))]),stringsAsFactors =FALSE)resumen[resumen$Train ==0,]```Posteriormente, eliminamos del conjunto test los casos, para los que no se puede crear un predictor```{r}test <- test[!(test$prov %in%c(5, 19)), ]test0 <- test[rowSums(is.na(test))==0,]```**Predicciones en el conjunto de test (modelos sin regularización)**```{r}#| code-fold: true#| code-summary: "Mostrar código"# Prediccionesp.Po.log <-predict(Po.log, test0 , type="response")p.Po.id <-predict(Po.id, test0 , type="response")p.BN.log <-predict(BN.log, test0 , type="response")p.ZI.Po <-predict(ZI.Po, test0 , type="response")p.ZI.BN <-predict(ZI.BN, test0 , type="response")nd <-data.frame(test0, of_var =log(test0$duracion_poliza/365.25)) p.Po.log.off <-predict(Po.log.off2, nd , type ="response")p.BN.log.off <-predict(BN.log.off, nd, type ="response")p.ZI.Po.off <-predict(ZI.Po.off, test0 , type="response")p.ZI.BN.off <-predict(ZI.BN.off, test0 , type="response")```Pasamos a generar las predicciones en el conjunto de prueba para los modelos sin regularización. Descartamos el modelo Poisson con pesos debido a que, como hemos podido observar, tiene un AIC significativamente más alto respecto al resto de modelos.```{r}#| code-fold: true#| code-summary: "Mostrar código"colores3 <-carto_pal(11, "RedOr")resultados <-data.frame(reales = test0$numerosi,Po_log = p.Po.log,Po_id = p.Po.id,BN_log = p.BN.log,ZI_Po = p.ZI.Po,ZI_BN = p.ZI.BN,Po_log_off = p.Po.log.off,BN_log_off = p.BN.log.off,ZI_Po_off = p.ZI.Po.off,ZI_BN_off = p.ZI.BN.off)resultados <-melt(resultados, id.vars ="reales", variable.name ="Modelo", value.name ="Predicciones")ggplot(resultados, aes(x = Predicciones, fill = Modelo)) +geom_density(alpha =0.4) +geom_density(data = resultados, aes(x = reales), fill ="black", alpha =0.2) +facet_wrap(~ Modelo) +scale_x_continuous(limits =c(0, 5)) +scale_fill_manual(values = colores3) +labs(x ="Número de Siniestros", y ="Densidad", title ="Comparación de Predicciones con Datos Reales") +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),axis.text.x =element_text(),panel.border =element_blank(),plot.margin =margin(2, 1, 1, 1, "cm"),legend.position ="none")```El gráfico muestra la comparación de las distribuciones de las predicciones obtenidas con distintos modelos en el test set, frente a los datos reales del mismo. Observamos que en la mayoría de modelos las predicciones siguen una tendencia parecida a la de los datos reales, con un pico inicial elevado cerca de 0 sinestros. Mientras que modelos como Po_id presentan diferencias marcadas respecto a realidad, otros como ZI_BN_off parecen ajustarse mejor.No obstante, para determinar correctamente cuál es el mejor modelo, procedemos a calcular los errores cuadráticos y absolutos.**Errores cuadráticos**```{r}#| code-fold: true#| code-summary: "Mostrar código"# Función de perdida cuadráticae.Po.log <-sum((p.Po.log - test0$numerosi)^2)e.Po.id <-sum((p.Po.id - test0$numerosi)^2)e.BN.log <-sum((p.BN.log - test0$numerosi)^2)e.ZI.Po <-sum((p.ZI.Po - test0$numerosi)^2)e.ZI.BN <-sum((p.ZI.BN - test0$numerosi)^2)e.Po.log.off <-sum((p.Po.log.off - test0$numerosi)^2)e.BN.log.off <-sum((p.BN.log.off - test0$numerosi)^2)e.ZI.Po.off <-sum((p.ZI.Po.off - test0$numerosi)^2)e.ZI.BN.off <-sum((p.ZI.BN.off - test0$numerosi)^2)e.base <-sum((mean(train$numerosi) - test0$numerosi)^2)``````{r}errores_cuadraticos <-c(e.Po.log, e.Po.id, e.BN.log, e.ZI.Po, e.ZI.BN, e.Po.log.off, e.BN.log.off, e.ZI.Po.off, e.ZI.BN.off, e.base)```El error cuadrático se utiliza para evaluar el ajuste de los modelos comparando las predicciones con los valores reales observados en el test set. Penaliza más fuertemente las grandes desviaciones entre las predicciones y los valores reales.**Errores valor absoluto**```{r}#| code-fold: true#| code-summary: "Mostrar código"a.Po.log <-sum(abs(p.Po.log - test0$numerosi))a.Po.id <-sum(abs(p.Po.id - test0$numerosi))a.BN.log <-sum(abs(p.BN.log - test0$numerosi))a.ZI.Po <-sum(abs(p.ZI.Po - test0$numerosi))a.ZI.BN <-sum(abs(p.ZI.BN - test0$numerosi))a.Po.log.off <-sum(abs(p.Po.log.off - test0$numerosi))a.BN.log.off <-sum(abs(p.BN.log.off - test0$numerosi))a.ZI.Po.off <-sum(abs(p.ZI.Po.off - test0$numerosi))a.ZI.BN.off <-sum(abs(p.ZI.BN.off - test0$numerosi))a.base <-sum(abs(mean(train$numerosi) - test0$numerosi))``````{r}errores_abs <-c(a.Po.log, a.Po.id, a.BN.log, a.ZI.Po, a.ZI.BN, a.Po.log.off, a.BN.log.off, a.ZI.Po.off, a.ZI.BN.off, a.base)```A diferencia del error cuadrático, el error absoluto no penaliza en exceso las grandes desviaciones, sino que simplemente calcula la suma de las diferencias absolutas entre las predicciones y los valores reales. Esto es útil cuando se quiere tener una idea más directa de la magnitud de los errores sin que los grandes errores afecten de forma desproporcionada el cálculo total.**Resumen errores (modelos sin regularización)**A continuación, se muestra una tabla resumiendo cuadráticos y absolutos obtenidos para cada modelo para facilitar la comparación entre los modelos ajustados.```{r}#| code-fold: true#| code-summary: "Mostrar código"errores_sin <-data.frame(modelo =c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", "Po.log.off", "BN.log.off", "ZI.Po.off", "ZI.BN.off", "Null model"), errores_cuadraticos, errores_abs)errores_sin[, 2:3] <-round(errores_sin[, 2:3], 2)errores_sin %>%kbl(align =rep('c', 3)) %>%kable_styling()``````{r}#| code-fold: true#| code-summary: "Mostrar código"errores_combinados <-rbind( errores_sin$errores_cuadraticos, errores_sin$errores_abs)colnames(errores_combinados) <-c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", "Po.log.off", "BN.log.off", "ZI.Po.off", "ZI.BN.off", "Null model")barplot(errores_combinados, beside =TRUE, names.arg =c("Po.log", "Po.id", "BN.log", "ZI.Po", "ZI.BN", "Po.log.off", "BN.log.off", "ZI.Po.off", "ZI.BN.off", "Null model"), col =c(colores2[2], colores2[7]),border ="white", main ="Comparación de Errores Cuadráticos y Absolutos", xlab ="Modelo", ylab ="Error", las =1, cex.names =0.6,cex.axis =0.8, cex.main =1.2)```Atendiendo a los errores y a las características de los modelos, elegiríamos el modelo de la Binomial Negativa. Como se ha comentado anteriormente, la distribución Binomial Negativa se ajusta mejor a nuestra muestra y es útil a la hora de sobrellevar la sobredispersión en los datos que se genera comúnmente en variables donde la varianza supera la media. A pesar de que los errores más bajos los presenta el modelo Poisson con enlace log, los valores tanto del error cuadrático como del absoluto son muy cercanos entre ambos modelos.Respecto a la evaluación clásica, aunque evaluar el AIC puede resultar en modelos más complejos con menos sobreajuste, los errores nos ofrecen una visión directa del rendimiento predictivo del modelo. Esto puede hacer que un modelo con un AIC más bajo pero con peores errores no será necesariamente la mejor opción si nos fijamos en la precisión de las predicciones.### 4.2 Modelos con regularizaciónTras una comparación inicial de modelos sin regularización basándonos en el ajuste de los datos, ahora incluiremos la regularización, penalizando la complejidad del modelo y promoviendo soluciones más simples y robustas.**Los Paquetes glmnet y mpath**Cuestiones sobre las funciones - Conversión de data.frame en matriz (construyendo variables dummy para factores e intercepto), necesario para predecir.model.matrix(\~., data)- Modelo con respuesta Poisson cv.glmnet(x, y, family="poisson", alpha, nfolds) glmnet(x, y, family="poisson", alpha) cvglmreg(formula, data, family="poisson", penalty, nfolds) glmreg(formula, data, family="poisson", penalty)- Modelo con respuesta Binomial Negativa cvglmregNB(formula, data, penalty, nfolds) glmregNB(formula, data, penalty)- Modelos inflados de ceros cvzipath(formula, data, family, link, penalty, nfolds) zipath(formula, data, family, link, penalty)- Predicción predict(object, newdata, type) predict(object, newx, which, type)Más información https://web.stanford.edu/\~hastie/glmnet/glmnet_alpha.html#poi https://cran.r-project.org/web/packages/mpath/mpath.pdfConversión a objetos matriz de predictores```{r}# Convertimos los data.frame en objetos matriz respetando factorestrain.m <-model.matrix(~., data=train)train.m <- train.m[,colnames(train.m) !="numerosi"]test.m <-model.matrix(~., data=test0)test.m <- test.m[,colnames(test.m) !="numerosi"]```**Modelos con regularización**Se estima un modelo, existe una solución, por cada valor del parámetro de penalización.::: {.callout-warning title="Descarte de modelos"}A lo largo del resto del trabajo, hemos decidido descartar algunos modelos debido a errores técnicos que no se pudieron solucionar durante el proceso. :::#### 4.2.1 Ajustes con glmnet```{r}#| code-fold: true#| code-summary: "Mostrar código"# Respuesta NormalNo.re <-glmnet(x = train.m, y = train$numerosi, alpha=0.5)# Respuesta PoissonPo.re <-glmnet(x = train.m, y = train$numerosi,family ="poisson", alpha =0.5)# Respuesta Poisson con offsetPo.re.off <-glmnet(x = train.m, y = train$numerosi, alpha =0.5, offset =log(train$duracion_poliza/365.25), family ="poisson")```Utilizamos glmnet para ajustar modelos con regularización mediante la penalización Elastic Net (alpha = 0.5), que combina las propiedades de las penalizaciones Lasso y Ridge.#### 4.2.2 Ajustes con mpath```{r}#| code-fold: true#| code-summary: "Mostrar código"# Respuesta PoissonPo.re2 <-glmreg(numerosi ~ . , data = train,family ="poisson", penalty="enet")# Respuesta Binomial NegativaBN.re <-glmregNB(numerosi ~ . , data = train, penalty="enet")# Respuestas Infladas de CerosZI.Po.re <-zipath(numerosi ~ .| . , data= train,family ="poisson", link ="logit", penalty="enet")ZI.BN.re <-zipath(numerosi ~ .| . , data= train,family ="negbin", link ="logit", penalty="enet")# Modelos Po y BN con offsetformula3 <- numerosi ~ potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia + g.anticia + edad2 + edad3# Modelos con offsetPo.re2.off <-glmreg(formula3, data = train, family ="poisson",offset =log(duracion_poliza/365.25), penalty="enet")BN.re.off <-glmregNB(formula3, data = train, offset =log(duracion_poliza/365.25), penalty="enet")# Modelos ZI con offsetformula4 <-formula(numerosi ~ potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia + g.anticia + edad2 + edad3 | potencia + antivehi + tipovehi + plazasvehi + edad + sexocondu + anticia + prov + usovehi2 + g.potencia + g.anticia + edad2 + edad3)ZI.Po.re.off <-zipath(formula4, data= train, family ="poisson", offset =log(duracion_poliza/365.25), link ="logit", penalty="enet", maxit =100)ZI.BN.re.off <-zipath(formula4, data= train, family ="negbin", offset =log(duracion_poliza/365.25), link ="logit", penalty="enet", maxit =100)```Ajustamos los modelos con mpath con fórmulas específicas para asegurarnos de que las variables predictoras relevantes estén correctamente incluidas en los ajustes. Este enfoque nos permite comparar diferentes técnicas y distribuciones, considerando características particulares de los datos y explorando la influencia de la regularización en cada caso.#### 4.2.3 Elección del modeloDisponemos de un modelo para cada lambda testado, ¿cuál elegimos?Vamos a considerar dos alternativas:- Predecimos con el modelo con el lambda que genera el "mejor" ajuste con medida clásica (AIC). - Elegimos el mejor modelo por Cross-Validation **Elección de modelos con regularización basándonos en el AIC y predicción**```{r}#| code-fold: true#| code-summary: "Mostrar código"# Modelo de respuesta Normalselec <-which(No.re$dev.ratio ==max(No.re$dev.ratio))p.No.re <-predict(No.re, test.m, type ="response")[, selec]# Modelo de respuesta Poisson (glmnet)selec <-which.max(Po.re$dev.ratio)p.Po.re <-predict(Po.re, test.m, type ="response")[, selec]# Modelo de respuesta Poisson y offset glmnetselec <-which.max(Po.re.off$dev.ratio)p.Po.re.off <-predict(Po.re.off, test.m, newoffset =log(test$duracion_poliza /365.25),type ="response")[, selec]# Modelo de respuesta Poisson (mpath)selec <-which.min(Po.re2$resdev)p.Po.re2 <-predict(Po.re2, test0, which = selec, type ="response")# Modelo de respuesta Poisson y offset (mpath)selec <-which.min(Po.re2.off$resdev)p.Po.re2.off <-exp(predict(Po.re2.off, test0, which = selec,newoffset =log(test0$duracion_poliza /365.25)))# Modelo de respuesta BN (mpath)selec <-which.min(AIC(BN.re))p.BN.re <-predict(BN.re, test0, which = selec, type ="response")# Modelo de respuesta BN y offset (mpath)selec <-which.min(AIC(BN.re.off))p.BN.re.off <-exp(predict(BN.re.off, test0, which = selec,newoffset =log(test0$duracion_poliza /365.25)))# Modelo de respuesta ZI-Po (mpath)selec <-which.min(AIC(ZI.Po.re))p.ZI.Po.re <-predict(ZI.Po.re, test0, which = selec,type ="response")# Modelo de respuesta ZI-Po y offset (mpath)selec <-which.min(AIC(ZI.Po.re.off))p.ZI.Po.re.off <-predict(ZI.Po.re.off, test0, which = selec,type ="response")# Modelo de respuesta ZI-BN (mpath)selec <-which.min(AIC(ZI.BN.re))p.ZI.BN.re <-predict(ZI.BN.re, test0, which = selec, type ="response")# Modelo de respuesta ZI-BN y offset (mpath)selec <-which.min(AIC(ZI.BN.re.off))p.ZI.BN.re.off <-predict(ZI.BN.re.off, test0, which = selec, type ="response")```Realizamos la selección y predicción de varios modelos ajustados con regularización, considerando qué tanto del comportamiento de los datos explica el modelo, o qué tan bueno es el modelo según el criterio del AIC. Además, incorporamos ajustes con offset para manejar correctamente la exposición al riesgo, asegurando que los modelos sean más interpretables y relevantes para el contexto.```{r}#| code-fold: true#| code-summary: "Mostrar código"resultados <-data.frame(reales = test0$numerosi,No_re = p.No.re,Po_re = p.Po.re,Po_re_off = p.Po.re.off,Po_re2 = p.Po.re2,BN_re = p.BN.re,ZI_Po_re = p.ZI.Po.re,Po_re2_off = p.Po.re2.off,BN_re_off = p.BN.re.off,ZI_Po_re_off = p.ZI.Po.re.off,ZI_BN_re = p.ZI.BN.re,ZI_BN_re_off = p.ZI.BN.re.off)resultados <-melt(resultados, id.vars ="reales", variable.name ="Modelo", value.name ="Predicciones")ggplot(resultados, aes(x = Predicciones, fill = Modelo)) +geom_density(alpha =0.4) +geom_density(data = resultados, aes(x = reales), fill ="black", alpha =0.2) +facet_wrap(~ Modelo) +scale_x_continuous(limits =c(0, 5)) +scale_fill_manual(values = colores3) +labs(x ="Número de Siniestros", y ="Densidad", title ="Comparación de Predicciones con Datos Reales") +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),axis.text.x =element_text(),panel.border =element_blank(),plot.margin =margin(2, 1, 1, 1, "cm"),legend.position ="none")```Comparamos las densidades de los predicciones de varios modelos con las de los datos reales del número de siniestros. Los modelos BN_re y BN_re_off destacan por ajustarse mejor a los datos reales, especialmente en los valores más frecuentes (0 y 1 siniestros). En cambio, modelos como Po_re_off tienden a sobreestimar los siniestros altos, lo que indica un ajuste menos preciso. Por otro lado, el modelo No_re, aunque tiene bajos errores, subestima frecuencias importantes y no representa bien la variabilidad de los datos. En general, los modelos basados en la Binomial Negativa logran capturar mejor las características reales de la distribución.Ahora, al igual que hicimos antes con los modelos sin regularización, comprobamos los errores cuadráticos y absolutos.**Errores cuadráticos**```{r}#| code-fold: true#| code-summary: "Mostrar código"# Funcion de perdida cuadráticae.No.re <-sum((p.No.re - test0$numerosi)^2)e.Po.re <-sum((p.Po.re - test0$numerosi)^2)e.Po.re.off <-sum((p.Po.re.off - test0$numerosi)^2)e.Po.re2 <-sum((p.Po.re2 - test0$numerosi)^2)e.BN.re <-sum((p.BN.re - test0$numerosi)^2)e.ZI.Po.re <-sum((p.ZI.Po.re - test0$numerosi)^2)e.ZI.BN.re <-sum((p.ZI.BN.re - test0$numerosi)^2)e.Po.re2.off <-sum((p.Po.re2.off - test0$numerosi)^2)e.BN.re.off <-sum((p.BN.re.off - test0$numerosi)^2)e.ZI.Po.re.off <-sum((p.ZI.Po.re.off - test0$numerosi)^2)e.ZI.BN.re.off <-sum((p.ZI.BN.re.off - test0$numerosi)^2)``````{r}errores_cuadraticos_re_AIC <-c(e.No.re, e.Po.re, e.Po.re.off, e.Po.re2, e.BN.re, e.ZI.Po.re, e.ZI.BN.re, e.Po.re2.off, e.BN.re.off, e.ZI.Po.re.off, e.ZI.BN.re.off)```**Errores valor absoluto**```{r}#| code-fold: true#| code-summary: "Mostrar código"# Funcion de perdida valor absolutoa.No.re <-sum(abs(p.No.re - test0$numerosi))a.Po.re <-sum(abs(p.Po.re - test0$numerosi))a.Po.re.off <-sum(abs(p.Po.re.off - test0$numerosi))a.Po.re2 <-sum(abs(p.Po.re2 - test0$numerosi))a.BN.re <-sum(abs(p.BN.re - test0$numerosi))a.ZI.Po.re <-sum(abs(p.ZI.Po.re - test0$numerosi))a.ZI.BN.re <-sum(abs(p.ZI.BN.re - test0$numerosi))a.Po.re2.off <-sum(abs(p.Po.re2.off - test0$numerosi))a.BN.re.off <-sum(abs(p.BN.re.off - test0$numerosi))a.ZI.Po.re.off <-sum(abs(p.ZI.Po.re.off - test0$numerosi))a.ZI.BN.re.off <-sum(abs(p.ZI.BN.re.off - test0$numerosi))``````{r}errores_abs_re_AIC <-c(a.No.re, a.Po.re, a.Po.re.off, a.ZI.BN.re, a.Po.re2, a.BN.re, a.ZI.Po.re, a.Po.re2.off, a.BN.re.off, a.ZI.Po.re.off, a.ZI.BN.re.off)```**Resumen errores (modelos con regularización) y AIC**```{r}#| code-fold: true#| code-summary: "Mostrar código"errores_re_AIC <-data.frame(modelo =c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", "BN.re", "ZI.Po.re", "Po.re2.off","BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off"), errores_cuadraticos_re_AIC, errores_abs_re_AIC)errores_re_AIC[, 2:3] <-round(errores_re_AIC[, 2:3], 2)errores_re_AIC %>%kbl(align =rep('c', 3)) %>%kable_styling()``````{r}#| code-fold: true#| code-summary: "Mostrar código"errores_combinados <-rbind( errores_re_AIC$errores_cuadraticos_re_AIC, errores_re_AIC$errores_abs_re_AIC)colnames(errores_combinados) <-c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", "BN.re", "ZI.Po.re", "Po.re2.off","BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off")barplot(errores_combinados, beside =TRUE, names.arg =c("No.re", "Po.re", "Po.re.off", "ZI.BN.re", "Po.re2", "BN.re", "ZI.Po.re", "Po.re2.off","BN.re.off", "ZI.Po.re.off", "ZI.BN.re.off"), col =c(colores2[2], colores2[7]),border ="white", main ="Comparación de Errores Cuadráticos y Absolutos", xlab ="Modelo", ylab ="Error", las =1, cex.names =0.6,cex.axis =0.8, cex.main =1.2)```Elegiríamos el modelo de la Binomial Negativa con offset. Entre los modelos evaluados, este modelo es el que presenta los mejores resultados generales, con el mínimo error tanto absoluto (163.79), como cuadrático (213.78). El valor del error cuadrático indica que el modelo maneja eficazmente posibles desviaciones extremas en los datos, mientras que el error absoluto refleja que en promedio sus predicciones se acercan más a los valores reales comparado con el resto de modelos.Comparando con los modelos sin regularización, observamos que los errores de los modelos con regularización son más bajos en general. Existe una disminución de los errores en el modelo regularizado que hemos elegido respecto al modelo sin regularización anteriormente preferido.#### 4.2.4 Elección de modelos usando cross-validation```{r}#| code-fold: true#| code-summary: "Mostrar código"train.m <-model.matrix(~., data=train)train.m <- train.m[,colnames(train.m) !="numerosi"]test.m <-model.matrix(~., data=test0)test.m <- test.m[,colnames(test.m) !="numerosi"]# Convertimos a matriz dispersatrain.m <-as(train.m, "dgCMatrix")# Modelo de respuesta Normalcv.No <-cv.glmnet(x=train.m, y=train$numerosi, alpha=0.5)# Modelo de respuesta Poisson (glmnet)cv.Po <-cv.glmnet(x=train.m, y=train$numerosi, alpha=0.5, family ="poisson")# Modelo de respuesta Poisson y offset (glmnet)cv.Po.off <-cv.glmnet(x = train.m, y = train$numerosi, alpha =0.5, family ="poisson", offset =log(train$duracion_poliza/365.25))# Modelo de respuesta Poisson (mpath)#cv.Po2 <- cv.glmreg(numerosi ~ . , data = train, family="poisson",# penalty = "enet")# Modelo de respuesta Poisson y offset (mpath)# Ajustar el modelo Poisson con penalización#cv.Po2.off <- cv.glmreg(formula3 , data = train_log, # family = "poisson", penalty="enet",# offset = log(duracion_poliza/365.25))# Modelo de respuesta BN (mpath)#cv.BN <- cv.glmregNB(formula3 , data= train, # penalty="enet", nfolds = 3, trace = TRUE)# Modelo de respuesta BN y offset (mpath)#cv.BN.off <- cv.glmregNB(formula3, data = train, # offset = log(train$duracion_poliza/365.25))# Modelo de respuesta ZI-Po (mpath)#cv.ZI.Po <- cv.zipath(numerosi ~ .| . , data= train,# family="poisson", link = "logit",# penalty="enet")# Modelo de respuesta ZI-BN (mpath)#cv.ZI.BN <- cv.zipath(numerosi ~ .| . , data= train, # family="negbin", link = "logit",# penalty="enet")# Modelo de respuesta ZI-Po y offset (mpath)# Variables y objetos auxiliarestrain$duracion_poliza.a <- train$duracion_poliza/365.25test0$duracion_poliza.a <- test0$duracion_poliza/365.25formula5 <-formula(numerosi~ potencia + antivehi+ tipovehi+ plazasvehi+ edad + sexocondu+ anticia+ prov+ usovehi2 + g.potencia+ g.anticia+ edad2 + edad3 + edad_discretizada+offset(log(duracion_poliza.a))| potencia + antivehi+ tipovehi+ plazasvehi+ edad + sexocondu+ anticia+ prov+ usovehi2 + g.potencia+ g.anticia+ edad2 + edad3 + edad_discretizada+offset(log(duracion_poliza.a)))#cv.ZI.Po.off<-cv.zipath(formula5, data= train, link = "logit",# family= "poisson", penalty="enet")# Modelo de respuesta ZI-BN y offset (mpath)#cv.BN.Po.off <- cv.zipath(formula5, data= train, family= "negbin",# link = "logit", penalty="enet")```El enfoque de cross-validation nos permite comparar el rendimiento de los modelos en diferentes conjuntos de entrenamiento y prueba, ayudando a identificar el modelo más adecuado según las métricas de error generadas durante la validación.**Predicciones con modelos ML con selección por CV**```{r}#| code-fold: true#| code-summary: "Mostrar código"p.cv.No <-predict(cv.No, test.m)p.cv.Po <-predict(cv.Po, test.m, type="response")p.cv.Po.off <-predict(cv.Po.off, test.m, type="response", newoffset =log(test$duracion_poliza/365.25))#p.cv.Po2 <-predict(cv.Po2$fit, # newx = test.m[, colnames(test.m) != "(Intercept)"], # which= cv.Po2$lambda.which, type = "response")#p.cv.Po2 <-predict(cv.Po2$fit, newx = test0)#p.cv.Po2.off <- predict(cv.Po2$fit, # newx = test.m[, colnames(test.m) != "(Intercept)"], # which= cv.Po2$lambda.which, type = "response", # newoffset = log(test0$duracion_poliza/365.25))#p.cv.BN <-exp(predict(cv.BN, test0))#p.cv.BN.off <- exp(predict(cv.BN.off, newx = test0, # newoffset = log(test0$duracion_poliza/365.25)))#p.cv.ZI.Po<-predict(cv.ZI.Po, test0, # which= cv.ZI.Po$lambda.which, # type="response")#p.cv.ZI.BN <-predict(ZI.BN.re, test0, # which= cv.ZI.BN$lambda.which, # type="response")p.cv.ZI.Po.off<-predict(ZI.Po.off, newdata = test0, which=cv.ZI.Po.off$lambda.which,type="response")p.cv.ZI.BN.off<-predict(ZI.BN.off, newdata = test0, which=cv.ZI.BN.off$lambda.which,type="response")```Tras entrenar los modelos utilizando técnicas de selección de hiperparámetros mediante validación cruzada, generamos las predicciones de los mismos.```{r}#| code-fold: true#| code-summary: "Mostrar código"resultados <-data.frame(reales = test0$numerosi,cv_No = p.cv.No,cv_Po = p.cv.Po,cv_Po_off = p.cv.Po.off,cv_ZI_Po_off = p.cv.ZI.Po.off,cv_ZI_BN_off = p.cv.ZI.BN.off)resultados <-melt(resultados, id.vars ="reales", variable.name ="Modelo", value.name ="Predicciones")resultados$Modelo <-recode(resultados$Modelo,"lambda.1se"="cv_No","lambda.1se.1"="cv_Po","lambda.1se.2"="cv_Po_off")ggplot(resultados, aes(x = Predicciones, fill = Modelo)) +geom_density(alpha =0.4) +geom_density(data = resultados, aes(x = reales), fill ="black", alpha =0.2) +facet_wrap(~ Modelo) +scale_x_continuous(limits =c(0, 5)) +scale_fill_manual(values = colores3) +labs(x ="Número de Siniestros", y ="Densidad", title ="Comparación de Predicciones con Datos Reales") +theme_minimal() +theme(text =element_text(size =12),plot.title =element_text(hjust =0.5, size =14, face ="bold"),axis.title.x =element_text(size =12), axis.title.y =element_text(size =12),axis.text.x =element_text(),panel.border =element_blank(),plot.margin =margin(2, 1, 1, 1, "cm"),legend.position ="none")```Al igual que antes, el gráfico compara la densidad de las predicciones, esta vez de los modelos validados con validación cruzada. El modelo Binomial con offset muestra el mejor ajuste, logrando capturar de forma precisa la distribución de los datos reales, especialmente en los siniestros más frecuentes. En contraste, cv_No subestima los valores intermedios, y cv_Po_off tiende a sobreestimar siniestros altos, reduciendo su precisión. Los resultados refuerzan la preferencia por cv_ZI_BN_off, que demuestra una mejor capacidad de generalización.**Errores cuadráticos**```{r}#| code-fold: true#| code-summary: "Mostrar código"e.cv.No <-sum(abs(p.cv.No - test0$numerosi)^2)e.cv.Po <-sum(abs(p.cv.Po - test0$numerosi)^2)e.cv.Po.off <-sum(abs(p.cv.Po.off - test0$numerosi)^2)#e.cv.Po2 <- sum(abs(p.cv.Po2 - test0$numerosi)^2)#e.cv.Po2.off <- sum(abs(p.cv.Po2.off - test0$numerosi)^2)#e.cv.BN <- sum(abs(p.cv.BN - test0$numerosi)^2)#e.cv.BN.off <- sum(abs(p.cv.BN.off - test0$numerosi)^2)#e.cv.ZI.Po <- sum(abs(p.cv.ZI.Po - test0$numerosi)^2)e.cv.ZI.Po.off <-sum(abs(p.cv.ZI.Po.off - test0$numerosi)^2)#e.cv.ZI.BN <- sum(abs(p.cv.ZI.BN - test0$numerosi)^2)e.cv.ZI.BN.off <-sum(abs(p.cv.ZI.BN.off - test0$numerosi)^2)``````{r}errores_cuadraticos_re_CV <-c(e.cv.No, e.cv.Po, e.cv.Po.off, e.cv.ZI.Po.off, e.cv.ZI.BN.off)```**Errores valor absoluto**```{r}#| code-fold: true#| code-summary: "Mostrar código"# Funcion de perdida valor absolutoa.cv.No <-sum(abs(p.cv.No - test0$numerosi))a.cv.Po <-sum(abs(p.cv.Po - test0$numerosi))a.cv.Po.off <-sum(abs(p.cv.Po.off - test0$numerosi))#a.cv.Po2 <- sum(abs(p.cv.Po2 - test0$numerosi))#a.cv.Po2.off <- sum(abs(p.cv.Po2.off - test0$numerosi))#a.cv.BN <- sum(abs(p.cv.BN - test0$numerosi))#a.cv.BN.off <- sum(abs(p.cv.BN.off - test0$numerosi))#a.cv.ZI.Po <- sum(abs(p.cv.ZI.Po - test0$numerosi))a.cv.ZI.Po.off <-sum(abs(p.cv.ZI.Po.off - test0$numerosi))#a.cv.ZI.BN <- sum(abs(p.cv.ZI.BN - test0$numerosi))a.cv.ZI.BN.off <-sum(abs(p.cv.ZI.BN.off - test0$numerosi))``````{r}errores_abs_re_CV <-c(a.cv.No, a.cv.Po, a.cv.Po.off, a.cv.ZI.Po.off, a.cv.ZI.BN.off)```**Resumen errores (modelos con regularización) y CV**```{r}#| code-fold: true#| code-summary: "Mostrar código"errores_re_CV <-data.frame(modelo =c("No.re", "Po.re", "Po.re.off","ZI.Po.re.off", "ZI.BN.re.off"), errores_cuadraticos_re_CV, errores_abs_re_CV)errores_re_CV[, 2:3] <-round(errores_re_CV[, 2:3], 2)errores_re_CV %>%kbl(align =rep('c', 3)) %>%kable_styling()``````{r}errores_combinados <-rbind( errores_re_CV$errores_cuadraticos_re_CV, errores_re_CV$errores_abs_re_CV)colnames(errores_combinados) <-c("No.re", "Po.re", "Po.re.off","ZI.Po.re.off", "ZI.BN.re.off")``````{r}#| code-fold: true#| code-summary: "Mostrar código"barplot(errores_combinados, beside =TRUE, names.arg =c("No.re", "Po.re", "Po.re.off","ZI.Po.re.off", "ZI.BN.re.off"), col =c(colores2[2], colores2[7]),border ="white", main ="Comparación de Errores Cuadráticos y Absolutos", xlab ="Modelo", ylab ="Error", las =1, cex.names =1,cex.axis =0.8, cex.main =1.2)```**Resumen errores modelos con regularización**```{r}#| code-fold: true#| code-summary: "Mostrar código"errores_re <-cbind(errores_re_AIC[c(1:3,9,8),], errores_re_CV[,-1])errores_re %>%kbl(align =rep('c', 3)) %>%kable_styling()```Volvemos a elegir el modelo BN.re.off como la mejor opción. Aunque fijándonos únicamente en el valor de los errores el modelo No.re podría ser más adecuado, este no incluye mecanismos para manejar sobredispersión ni incorpora términos de regularización que ayuden a evitar el sobreajusteEl modelo de la Binomial Negativa es particularmente adecuado para nuestro conjunto de datos, ya que, como se ha mencionado, maneja eficientemente la sobredispersión que caracteriza a la variable dependiente. Además, la regularización, a través del término offset, contribuye a evitar problemas de sobreajuste y mejora la estabilidad del modelo.Comparando con los modelos anteriores, BN.re.off demuestra una mejora significativa respecto a los modelos sin regularización, como BN.log.off, que aunque presenta buenos resultados, no alcanza el desempeño del modelo seleccionado en las métricas clave. En relación a los modelos evaluados únicamente con AIC (evaluación clásica), BN.re.off logra un balance superior entre ajuste y capacidad predictiva, superando a opciones como ZI.Po.off y ZI.BN, cuyos valores de AIC más bajos no se traducen en un mejor desempeño en validación cruzada.## 5. ConclusionesEn conclusión, dado que nuestro estudio se basa en una muestra donde la varianza de los datos supera la media, la Binomial Negativa es preferible gracias a su capacidad para capturar la estructura de estos datos, manejar la sobredispersión y garantizar predicciones confiables en nuevos conjuntos.La decisión de reducir la muestra nos permitió agilizar la convergencia de los modelos y continuar con el análisis. Sin embargo, es importante destacar que esto podría haber limitado la capacidad de los modelos para captar patrones más complejos presentes en la muestra completa. Adicionalmente, se eliminaron algunas observaciones específicas para evitar problemas al realizar las predicciones, como aquellas correspondientes a niveles de factores no presentes en el conjunto de entrenamiento. Estos ajustes fueron necesarios para garantizar la ejecución, pero podrían haber afectado la representatividad de los resultados. Por lo tanto, sería interesante repetir el análisis utilizando la totalidad de los datos originales en un equipo de soporte más preparado para manejar tal cantidad de datos.## 6. Referencias\[Callout Blocks – quarto. (n.d.). Quarto.\] (https://quarto.org/docs/authoring/callouts.html)\[INE - Instituto Nacional de Estadística. (n.d.). INEbase/ Clasificaciones /Relación de municipios, provincias, comunidades y ciudades autónomas con sus códigos / Relación de comunidades y ciudades autónomas con sus códigos. INE - Instituto Nacional De Estadística.\] (https://www.ine.es/daco/daco42/codmun/cod_ccaa_provincia.htm)\[¿Qué es la correlación estadística y cómo interpretarla? (n.d.). Máxima Formación.\] (https://www.maximaformacion.es/blog-dat/que-es-la-correlacion-estadistica-y-como-interpretarla/)\[colaboradores de Wikipedia. (2023, November 22). Distribución log-normal. Wikipedia, La Enciclopedia Libre.\] (https://es.wikipedia.org/wiki/Distribuci%C3%B3n_log-normal)\[Color palettes – Data Visualization with R. (n.d.). Data Visualization With R.\] (https://datavizs24.classes.andrewheiss.com/resource/colors.html)\[What are continuous probability distributions & their 8 common types? \| KNIME. (n.d.). KNIME.\] (https://www.knime.com/blog/continuous-probability-distribution)
